An improved quasi-reversibility method for a terminal-boundary value multi-species model with white Gaussian noise

TL;DR

This paper advances the quasi-reversibility method for stochastic terminal-boundary parabolic problems, providing convergence analysis, error bounds, and numerical validation for a multi-species reaction-diffusion model with white Gaussian noise.

Contribution

It introduces an improved quasi-reversibility approach with convergence analysis and error estimates for stochastic multi-species models, extending previous deterministic methods.

Findings

Established logarithmic stability estimates for the stabilized operator.

Derived error bounds for concentrations and gradients under stochastic noise.

Numerical examples confirm theoretical error estimates.

Abstract

Upon the recent development of the quasi-reversibility method for terminal value parabolic problems in \cite{Nguyen2019}, it is imperative to investigate the convergence analysis of this regularization method in the stochastic setting. In this paper, we positively unravel this open question by focusing on a coupled system of Dirichlet reaction-diffusion equations with additive white Gaussian noise on the terminal data. In this regard, the approximate problem is designed by adding the so-called perturbing operator to the original problem and by exploiting the Fourier reconstructed terminal data. By this way, Gevrey-type source conditions are included, while we successfully maintain the logarithmic stability estimate of the corresponding stabilized operator, which is necessary for the error analysis. As the main theme of this work, we prove the error bounds for the concentrations and for the concentration gradients, driven by a large amount of weighted energy-like controls involving the expectation operator. Compared to the classical error bounds in $L^2$ and $H^1$ that we obtained in the previous studies, our analysis here needs a higher smoothness of the true terminal data to ensure their reconstructions from the stochastic fashion. Two numerical examples are provided to corroborate the theoretical results.

Equations

$$\{\begin{array}{ll}{u}_t-{\mathcal{D}}_1\left({\ell }_0\left(u\right)\left(t\right)\right)\Delta u=F\left(x,t;u;v\right),& \\ v_t-{\mathcal{D}}_2\left({\ell }_0\left(v\right)\left(t\right)\right)\Delta v=G\left(x,t;u;v\right)& \text{for }\left(x,t\right)\in Q_T.\end{array}$$

$$u_f^{\epsilon }\left(x\right)=u_f\left(x\right)+\epsilon {\xi }_1\left(x\right),v_f^{\epsilon }\left(x\right)=v_f\left(x\right)+\epsilon {\xi }_2\left(x\right),$$

$${\mathcal{D}}_i\left({\ell }_0\left(w\right)\left(t\right)\right)={\mathcal{D}}_i\left({\int }_{\Omega }f(x)w(x,t)dx\right),{\ell }_0(w)(t):={\int }_{\Omega }f(x)w(x,t)dx,$$

$$u_f^{\epsilon ,j}$$

$$v_f^{\epsilon ,j}$$

$$0\le {\mu }_0<{\mu }_1\le {\mu }_2\le \dots ,\text{and }\underset{j\to \infty }{\lim }{\mu }_j=\infty .$$

$${\Vert {\mathbf{Q}}_{\epsilon }^{\beta }u\Vert }_{{\left[L^2\left(\Omega \right)\right]}^N}\le \frac{C_0}{\gamma \left(T,\beta \right)}{\Vert u\Vert }_{\mathbb{W}}\text{for any }u\in \mathbb{W}.$$

$${\Vert {\mathbf{P}}_{\epsilon }^{\beta }u\Vert }_{{\left[L^2\left(\Omega \right)\right]}^N}\le C_1\log \left(\gamma \left(T,\beta \right)\right){\Vert u\Vert }_{{\left[L^2\left(\Omega \right)\right]}^N}\text{for any }u\in {\left[L^2\left(\Omega \right)\right]}^N.$$

$$\gamma \left(T,\beta \right)\ge 1,\underset{\beta \to 0^{+}}{\lim }\gamma \left(t,\beta \right)=\infty \text{for all }t\in \left(0,T\right].$$

$$\underline{M}\le {\mathcal{D}}_i\left(\xi \right)\le M_1<\overline{M}\text{for any }\xi \in \mathbb{R}.$$

$$\left|F\left(x,t;u_1;v_1\right)-F\left(x,t;u_2;v_2\right)\right|+\left|G\left(x,t;u_1;v_1\right)-G\left(x,t;u_2;v_2\right)\right|$$

$$\le L\left(\ell \right)\left(\left|u_1-u_2\right|+\left|v_1-v_2\right|\right),$$

$$\left|{\mathcal{D}}_i\left({\ell }_0\left(u_1\right)\left(t\right)\right)-{\mathcal{D}}_i\left({\ell }_0\left(u_2\right)\left(t\right)\right)\right|\le \tilde{L}\Vert u_1\left(\cdot ,t\right)-u_2\left(\cdot ,t\right)\Vert \text{for }{u}_1,u_2\in L^2\left(\Omega \right).$$

$$L^2(\Omega )\ni U_f^{\epsilon ,n}\left(x\right):=\sum _{j=1}^n⟨u_f^{\epsilon },{\varphi }_j⟩{\varphi }_j\left(x\right),L^2(\Omega )\ni V_f^{\epsilon ,n}\left(x\right):=\sum _{j=1}^n⟨v_f^{\epsilon },{\varphi }_j⟩{\varphi }_j\left(x\right).$$

$$\mathbb{E}{\Vert U_f^{\epsilon ,n}-u_f\Vert }^2\le {\epsilon }^{2}n+\frac{{\Vert u_f\Vert }_{H^{2p}\left(\Omega \right)}^2}{{\mu }_n^{2p}},\mathbb{E}{\Vert V_f^{\epsilon ,n}-v_f\Vert }^2\le {\epsilon }^{2}n+\frac{{\Vert v_f\Vert }_{H^{2p}\left(\Omega \right)}^2}{{\mu }_n^{2p}},$$

$${\mathbf{Q}}_{\epsilon }^{\beta }u=\overline{M}\sum _{{\mu }_j>\frac{1}{\overline{M}T}\log \left(\gamma \left(T,\beta \right)\right)}{\mu }_j⟨u,{\varphi }_j⟩{\varphi }_j\text{for }u\in L^2\left(\Omega \right).$$

$${\mathbb{G}}_{p,q}=\left\{u\in L^2\left(\Omega \right):\sum _{j\in \mathbb{N}}{\mu }_j^q{e}^{2p{\mu }_j}{\left|⟨u,{\varphi }_j⟩\right|}^2<\infty \right\}.$$

$${⟨u_1,u_2⟩}_{{\mathbb{G}}_{p,q}}=⟨{\left(-\Delta \right)}^{q/2}{e}^{p\sqrt{-\Delta }}{u}_1,{\left(-\Delta \right)}^{q/2}{e}^{p\sqrt{-\Delta }}{u}_2⟩,$$

$${\Vert u_1\Vert }_{{\mathbb{G}}_{p,q}}=\sqrt{\sum _{j\in \mathbb{N}}{\mu }_j^q{e}^{2p{\mu }_j}{\left|⟨u,{\varphi }_j⟩\right|}^2}.$$

$${\mathbf{P}}_{\epsilon }^{\beta }u=-\overline{M}\sum _{{\mu }_j\le \frac{1}{\overline{M}T}\log \left(\gamma \left(T,\beta \right)\right)}{\mu }_j⟨u,{\varphi }_j⟩{\varphi }_j,$$

$$F_{{\ell }^{\epsilon }}\left(x,t;u,v\right):=\{\begin{array}{ll}F\left(x,t;{\ell }^{\epsilon };{\ell }^{\epsilon }\right)& \text{if }\max \left\{u,v\right\}>{\ell }^{\epsilon },\\ F\left(x,t;u;v\right)& \text{if }\max \left\{u,v\right\}\in \left[-{\ell }^{\epsilon },{\ell }^{\epsilon }\right],\\ F\left(x,t;-{\ell }^{\epsilon };-{\ell }^{\epsilon }\right)& \text{if }\max \left\{u,v\right\}<-{\ell }^{\epsilon },\end{array}$$

$$\left|F_{{\ell }^{\epsilon }}\left(x,t;u_1;v_1\right)-F_{{\ell }^{\epsilon }}\left(x,t;u_2;v_2\right)\right|+\left|G_{{\ell }^{\epsilon }}\left(x,t;u_1;v_1\right)-G_{{\ell }^{\epsilon }}\left(x,t;u_2;v_2\right)\right|$$

$$\le 2L\left({\ell }^{\epsilon }\right)\left(\left|u_1-u_2\right|+\left|v_1-v_2\right|\right),$$

$$\{\begin{array}{ll}{u}_t^{\epsilon ,n}-{\mathcal{D}}_1\left(u^{\epsilon ,n}\right)\Delta u^{\epsilon ,n}-{\mathbf{Q}}_{\epsilon }^{\beta }{u}^{\epsilon ,n}=F\left(x,t;u^{\epsilon ,n};v^{\epsilon ,n}\right),& \\ v_t^{\epsilon ,n}-{\mathcal{D}}_2\left(v^{\epsilon ,n}\right)\Delta v^{\epsilon ,n}-{\mathbf{Q}}_{\epsilon }^{\beta }{v}^{\epsilon ,n}=G\left(x,t;u^{\epsilon ,n};v^{\epsilon ,n}\right)& \text{for }\left(x,t\right)\in Q_T.\end{array}$$

$$\{\begin{array}{ll}{u}_t^{\epsilon ,n}+{\overline{\mathcal{D}}}_1\left(u^{\epsilon ,n}\right)\Delta u^{\epsilon ,n}=F_{{\ell }^{\epsilon }}\left(x,t;u^{\epsilon ,n};v^{\epsilon ,n}\right)+{\mathbf{P}}_{\epsilon }^{\beta }{u}^{\epsilon ,n},& \\ v_t^{\epsilon ,n}+{\overline{\mathcal{D}}}_2\left(v^{\epsilon ,n}\right)\Delta v^{\epsilon ,n}=G_{{\ell }^{\epsilon }}\left(x,t;u^{\epsilon ,n};v^{\epsilon ,n}\right)+{\mathbf{P}}_{\epsilon }^{\beta }{v}^{\epsilon ,n}& \text{for }\left(x,t\right)\in Q_T,\end{array}$$

$$u^{\epsilon ,n}\left(x,T\right)=U_f^{\epsilon ,n}\left(x\right),v^{\epsilon ,n}\left(x,T\right)=V_f^{\epsilon ,n}\left(x\right).$$

$$\{\begin{array}{l}{U}_t^{\epsilon ,n}+{\overline{\mathcal{D}}}_1\left(e^{{\rho }_{\epsilon }\left(T-t\right)}{U}^{\epsilon ,n}\right)\Delta U^{\epsilon ,n}-{\rho }_{\epsilon }{U}^{\epsilon ,n}\\ =e^{{\rho }_{\epsilon }\left(t-T\right)}{F}_{{\ell }^{\epsilon }}\left(e^{{\rho }_{\epsilon }\left(T-t\right)}{U}^{\epsilon ,n};e^{{\rho }_{\epsilon }\left(T-t\right)}{V}^{\epsilon ,n}\right)+{\mathbf{P}}_{\epsilon }^{\beta }{U}^{\epsilon ,n},\\ V_t^{\epsilon ,n}+{\overline{\mathcal{D}}}_2\left(e^{{\rho }_{\epsilon }\left(T-t\right)}{V}^{\epsilon ,n}\right)\Delta V^{\epsilon ,n}-{\rho }_{\epsilon }{V}^{\epsilon ,n}\\ =e^{{\rho }_{\epsilon }\left(t-T\right)}{G}_{{\ell }^{\epsilon }}\left(e^{{\rho }_{\epsilon }\left(T-t\right)}{U}^{\epsilon ,n};e^{{\rho }_{\epsilon }\left(T-t\right)}{V}^{\epsilon ,n}\right)+{\mathbf{P}}_{\epsilon }^{\beta }{V}^{\epsilon ,n},\end{array}$$

$$\left(U^{\epsilon ,n},V^{\epsilon ,n}\right)\in {\left[L^2\left(0,T;H_0^1\left(\Omega \right)\right)\cap L^{\infty }\left(0,T;L^2\left(\Omega \right)\right)\right]}^2$$

$$\frac{d}{dt}⟨U^{\epsilon ,n},{\psi }_1⟩$$

$$=e^{{\rho }_{\epsilon }\left(t-T\right)}⟨F_{{\ell }^{\epsilon }}\left(e^{{\rho }_{\epsilon }\left(T-t\right)}{U}^{\epsilon ,n};e^{{\rho }_{\epsilon }\left(T-t\right)}{V}^{\epsilon ,n}\right),{\psi }_1⟩+⟨{\mathbf{P}}_{\epsilon }^{\beta }{U}^{\epsilon ,n},{\psi }_1⟩,$$

$$\frac{d}{dt}⟨V^{\epsilon ,n},{\psi }_2⟩$$

$$=e^{{\rho }_{\epsilon }\left(t-T\right)}⟨G_{{\ell }^{\epsilon }}\left(e^{{\rho }_{\epsilon }\left(T-t\right)}{U}^{\epsilon ,n};e^{{\rho }_{\epsilon }\left(T-t\right)}{V}^{\epsilon ,n}\right),{\psi }_2⟩+⟨{\mathbf{P}}_{\epsilon }^{\beta }{V}^{\epsilon ,n},{\psi }_2⟩,$$

$$\underset{\epsilon \to 0^{+}}{\lim }{\gamma }^{C_{1}T}(T,\beta )\epsilon \sqrt{n}=K_0\in (0,\infty ),\underset{\epsilon \to 0^{+}}{\lim }{\gamma }^{C_{1}T}(T,\beta ){\mu }_n^{-p}=K_1\in (0,\infty ).$$

$$u,v\in C([0,T];L^2(\Omega ))\cap L^2(0,T;\bar{\mathbb{W}})\cap L^{\infty }(0,T;H_0^1(\Omega )\cap L^{\infty }(\Omega )),$$

Full text

{AMS}

62P10, 65J05, 65J20, 35K92, 60H35

1 Introduction

The main purpose of this paper is to carry out the error analysis of the recently proposed quasi-reversibility (QR) method in a stochastic setting for a class of terminal-boundary value multi-species model. As concluded in [26], the error bounds we obtained there can be very helpful in the finite element settings due to the variational framework we choose. Nonetheless, it still questions us about if this approach can be modified to somewhat get the convergence analysis in some specific stochastic setting involved in the partial differential equations (PDEs) we are chasing. In this work, we provide a positive answer by taking into account a nonlocal coupled system of nonlinear reaction-diffusion equations with additive white Gaussian noise on the terminal data. Besides, we significantly strengthen the applicability of this general QR framework by introducing a ready-to-use scheme from our computational standpoint.

1.1 Statement of the problem

Let $0<T<\infty$ be the final time of observation and $\Omega\subset\mathbb{R}^{d}$ for $\mathbb{N}^{*}\ni d\leq 3$ be the domain of interest, which is open, connected, and bounded with a sufficiently smooth boundary. Denoted by $Q_{T}=\Omega\times(0,T)$, this work is devoted to finding $u,v:\overline{Q_{T}}\to\mathbb{R}$ as two solutions of the following evolution system:

[TABLE]

To complete the terminal-boundary value problem we wish to solve, the system (1) is then supplemented with the terminal conditions $u(x,T)=u_{f}(x)$, $v(x,T)=v_{f}(x)$ and the Dirichlet boundary conditions $u=v=0$ on the boundary $\partial\Omega$. Suppose further that the terminal data are observed in the presence of white noise processes $\xi_{1}$, $\xi_{2}$ in terms of

[TABLE]

where $\varepsilon\in(0,1)$ is used for measuring the amplitude of the noise.

In this scenario, our attempt is to seek the initial values $u(x,0)=u_{0}(x)$ and $v(x,0)=v_{0}(x)$ when we only know the measurements of the terminal data $u_{f}$ and $v_{f}$. We accentuate that this seeking can only be done by regularization since our problem is, in general, exponentially unstable; cf. [4, 5, 11, 1, 20, 14] for a glimpse of statistical inverse problems. Starting from the source reconstruction model of diffusive, competitive and reactive brain tumor cells (cf. [13] and references cited therein), solving (1) is essential to understand the localization of the tumor source and then to possibly advance the treatment of brain tumors. This biological context conjures up the following images. Let $u$ and $v$ be the normal (healthy) and abnormally growing tissue cells densities in a brain region, respectively. In the perspective of (1), we assume that the movements of each kind of cells are dominantly influenced by the whole population of the corresponding type. Usually, this significant impact brings us to the nonlocal form of diffusion (see [7, 3, 9] for the mathematical background of this diffusion in the forward model), which reads as

[TABLE]

for $i=1,2$, $w\in\mathbb{R}$ and where $f$ is a given sufficiently smooth weight function. Largely motivated by the typical expression, we considered in [28] a slight generalization of this diffusion for the sole purpose of extending the applicability of the QR method we are studying.

The next concern that we would point out in the continuous model (1) lies in the nonlinearities $F$ and $G$ which can be considered as reaction, death and proliferation rates involved in the active network under scrutiny. In fact, we are able to suppose that at the low-grade gliomas regime, one knows that the brain tumor cells start growing slowly by infiltrating into the healthy brain cells and simultaneously, they misdirect the T cells of the immune system to avoid being beaten. Besides, we can assume that for some vigorous people, those gliomas are still identified at the beginning and then are very slightly removed by some reaction rates from the immune system. Accordingly, this entire incident is well-suited to the consideration of the coupled system (1), provided that the assumptions on these nonlinear terms will cover most of the real-world contexts. Needless to say, one may also detail several well-known chemical models contained in (1) such as the two-species Lotka–Volterra competition-diffusion model, the Brusselator coupled system for the autocatalytic reaction and those already mentioned in [26]. Eventually, we remark that even though the evolution system (1) can be expressed in a closed-form, it is better to contemplate the whole mathematical treatment for every single term in there.

1.2 Settings of the statistical terminal data

Cf. [4], we recall standard assumptions on the stochastic noise we want to address in the model. In the following, $\mathcal{H}$ is a Hilbert space.

Definition 1.1**.**

*The stochastic error is a Hilbert-space process, i.e. a bounded linear operator $\xi:\mathcal{H}\to L^{2}(\tilde{\Omega},\mathcal{A},P)$ where $(\tilde{\Omega},\mathcal{A},P)$ is a complete probability space and $L^{2}(\cdot)$ is the space of all square integrable measurable functions. *

This way we are able to define the random variables $\left\langle\xi,g_{j}\right\rangle$ for $j=1,2$ for all $g_{1},g_{2}\in\mathcal{H}$ by definition $\mathbb{E}\left\langle\xi,g_{j}\right\rangle=0$. Furthermore, we define its covariance $\text{Cov}_{\xi}$ as the bounded linear operator mapping from $\mathcal{H}$ onto itself such that $\left\langle\text{Cov}_{\xi}g_{1},g_{2}\right\rangle=\text{Cov}\left(\left\langle\xi,g_{1}\right\rangle,\left\langle\xi,g_{2}\right\rangle\right)$.

Definition 1.2**.**

*We say that $\xi$ is a white noise process in $\mathcal{H}$, if $\text{Cov}_{\xi}=I$ and the induced random variables are Gaussian: for all $g_{1},g_{2}\in\mathcal{H}$, the random variables $\left\langle\xi,g_{j}\right\rangle$ have distributions $\mathcal{N}\left(0,\left\|g_{j}\right\|^{2}\right)$ and $\text{Cov}\left(\left\langle\xi,g_{1}\right\rangle,\left\langle\xi,g_{2}\right\rangle\right)=\left\langle g_{1},g_{2}\right\rangle$. *

Accordingly, assume that the observations (2) can only be obtained in a discretized or binned form. This means that we only have vectors of normally distributed random variables $\left\{u_{f}^{\varepsilon,j}\right\}_{j=\overline{1,n}}$, $\left\{v_{f}^{\varepsilon,j}\right\}_{j=\overline{1,n}}$ given by

[TABLE]

where $n\in\mathbb{N}$ is the number of steps of discrete observations and $\phi_{j}$ is itself taken from the orthonormal basis $\left\{\phi_{j}\right\}_{j\in\mathbb{N}}$ of $L^{2}(\Omega)$. Recall that due to the smoothness of $\Omega$, the existence of this basis is guaranteed, where $\phi_{j}\in H_{0}^{1}(\Omega)\cap C^{\infty}(\overline{\Omega})$ solves the basic eigenvalue problem $-\Delta\phi_{j}(x)=\mu_{j}\phi_{j}(x)$ for $x\in\Omega$. Additionally, the Dirichlet eigenvalues $\left\{\mu_{j}\right\}_{j\in\mathbb{N}}$ form an infinite sequence which goes to infinity, viz.

[TABLE]

As a consequence, one can prove that $\left\langle\xi_{1},\phi_{j}\right\rangle$ and $\left\langle\xi_{2},\phi_{j}\right\rangle$ for $j=\overline{1,n}$ are i.i.d. standard Gaussian random variables.

1.3 QR-based methods and our novelty

The QR method has a long remarkable history since the pioneering monograph [21]. During the development of the QR method for inverse PDEs with deterministic noise, people usually focus on

• •

the spectral methods that allow us to consider the solution in a mild presentation (cf., e.g., [28, 8, 15]), but this typical method seems hard to handle the practical error control since one struggles with the Fourier accumulation for the nonlinear contexts;

• •

the Carleman-type estimate (cf., e.g., [19, 18]) that includes the uniqueness result and convergence analysis, but currently, it only works for the short time observation with the linear equation.

Cf. [26], our approach removes this fence by relying on the so-called perturbing and stabilized operators, benefited from the essential bound of the nonlinear diffusion coefficient.

Definition 1.3** (perturbing operator).**

The linear mapping $\mathbf{Q}_{\varepsilon}^{\beta}:[L^{2}(\Omega)]^{N}\to[L^{2}(\Omega)]^{N}$ is said to be a perturbing operator if there exist a function space $\mathbb{W}\subset[L^{2}(\Omega)]^{N}$ and an $\varepsilon$-independent constant $C_{0}>0$ such that

[TABLE]

Definition 1.4** (stabilized operator).**

The linear mapping $\mathbf{P}_{\varepsilon}^{\beta}:[L^{2}(\Omega)]^{N}\to[L^{2}(\Omega)]^{N}$ is said to be a stabilized operator if there exists an $\varepsilon$-independent constant $C_{1}>0$ such that

[TABLE]

In Definitions 1.3 and 1.4, $N$ denotes the number of species involved in the model and in this work, $N=2$. Meanwhile, we denote by $\beta=\beta(\varepsilon)\in(0,1)$ the regularization parameter satisfying $\lim_{\varepsilon\to 0^{+}}\beta(\varepsilon)=0$. The function $\gamma:[0,T]\times(0,1)\to\mathbb{R}$ indicates the decay behaviour of the perturbing operator, which stems from the source condition measuring the high smoothness of the true solution. In our analysis, we require: for any $\beta>0$, there holds

[TABLE]

In this study, we aim at developing this new QR approach and estimating its quality in a statistical setting. Cf. [4] the problem is well known to be exponentially ill-posed (compared to those mildly ill-posed introduced in the same reference); therefore, it is natural that we start from a very standard noise process. Yet, it is still a theoretical study that requires knowledge of the magnitude of the noise $\varepsilon$. Our theoretical analysis in this work and in all previous studies relies very much on the noise level as we want the (stably) approximate solution is close to the true one when $\varepsilon$ gets smaller. Without this analysis, there is nothing to ensure that the QR method performs well in solving such a highly challenging nonlinear inverse problem. Besides, this work would prepare a playground for the future evolution of our method. As we want to approach the real-world applications, it will concern more practical frameworks.

Although it is clear that in the framework of deterministic noise we can prove the error estimates for the QR method, it does not ensure that one can adapt those to the statistical inverse model in this work, including the references cited above. Besides, the existing literature on regularization of the system (1) is very limited due to the inception stage. Starting off with the linear version of (1) (where $\mathcal{D}_{i}$, $F$ and $G$ are independent of the solutions), we are aware of the presence of [23] where the authors exploited the trigonometric method in nonparametric regression to make use of the cut-off regularization in the statistical setting. This method is, however, very challenging in the nonlinear perspective as postulated in, e.g., [17]. It is worth noting in [17] that the authors essentially applied the QR-based method, where the Gevrey source condition is taken into account, to regularize a semilinear parabolic problem. Henceforth, our contribution herein can also be seen as an improvement of [17], since the nonlocal spatial operator becomes rather challenging from the numerical standpoint in [17].

From now on, some impedimenta that we will meet in proofs of our analysis should be revealed, except what we were very much concerned in [26]. Due to the white noise processes we are taking into account, our observations, i.e. the measured terminal data, do not belong to $L^{2}(\Omega)$, but acts on $L^{2}(\Omega)$ by virtue of (3)–(4). In principle, the random noise is large compared to the deterministic one. Thus, the adaptation of [26] to this statistical scenario is not straightforward. Aside from the high smoothness of the true solution, we need the true terminal function to be very smooth to obtain rates of convergence in expectation.

1.4 Outline of the paper

The rest of the paper is organized as follows. In Section 2, we introduce the notation and working assumptions for our analysis below. We also provide a “computable” example to validate the presence of the perturbation in Definition 1.3. In addition, we delineate a way to reconstruct the terminal data from the stochastic setting under consideration. Our main results are reported in Section 3. Based on the proposed QR framework, we establish a regularization for the time-reversed system (1); cf. Section 3.1. Convergence rates of the scheme are thoroughly explored in Section 3.2; cf. Theorems 3.3 and 3.4. Proofs of these error estimates are detailed in Sections 3.2.1 and 3.2.2, respectively. We close this main section by several discussions in Section 3.2.3. Finally, two numerical tests are provided in Section 4 to verify our convergence analysis, and the conclusions follow in Section 5.

2 Preliminaries

In the sequel, wherever $\left\langle\cdot,\cdot\right\rangle$ and $\left\|\cdot\right\|$ are present, we mean the $L^{2}$ inner product and its corresponding norm. Meanwhile, the other standard Sobolev norms will be specified, if used. Also, we use $\mathcal{D}_{i}(u)(t)$ as $\mathcal{D}_{i}(\ell_{0}(u)(t))$ for ease of presentation. In this statistical inverse problem, we assume that the true solution exists uniquely and its regularity is assumed in Theorems 3.3 and 3.4. The backward uniqueness result for this problem will be studied in the near future. To this end, we use the following assumptions:

$\left(\text{A}_{1}\right)$ For $i=1,2$, the measurable functions $\mathcal{D}_{i}>0$ is such that the mapping $\xi\mapsto\mathcal{D}_{i}\left(\xi\right)$ is continuous for $\xi\in\mathbb{R}$. Moreover, there exist $\varepsilon$-independent $\underline{M},M_{1},\overline{M}>0$ such that

[TABLE]

$\left(\text{A}_{2}\right)$ The source functions $F,G$ are measurable and locally Lipschitz-continuous in the sense that

[TABLE]

for $\max\left\{\left|u_{1}\right|,\left|u_{2}\right|,\left|v_{1}\right|,\left|v_{2}\right|\right\}\leq\ell$ for some $\ell>0$.

$\left(\text{A}_{3}\right)$ The true final conditions $u_{f},v_{f}$ belong to $H^{2p}(\Omega)$ for $p>0$.

$\left(\text{A}_{1}^{\prime}\right)$ For $i=1,2$, the diffusion $\mathcal{D}_{i}$ is globally Lipschitz-continuous in the sense that there exists $\varepsilon$-independent $\tilde{L}$ such that

[TABLE]

Since our observations $u_{f}^{\varepsilon}$ and $v_{f}^{\varepsilon}$ are not in general elements of $L^{2}(\Omega)$, we need the presence of $\left(\text{A}_{3}\right)$ to handle the following lemma (cf. [25]):

Lemma 2.1** (Truncated Fourier reconstruction).**

Given $n\in\mathbb{N}^{*}$ the number of steps of discrete observations, we define

[TABLE]

Then with the aid of $\left(\text{A}_{3}\right)$, one has

[TABLE]

*where $\mu_{n}>0$ is the $n$th Dirichlet eigenvalue of the Laplacian operator. *

The existence of the perturbing and stabilized operators (for Definitions 1.3 and 1.4) were commenced in [26], mimicking the stochastic gradient descent algorithm to obtain the function space $\mathbb{W}$ as a Gevrey111See again the Gevrey-like space defined in [26, Section 5.2]. class of real-analytic functions. Below, we exemplify another one using the truncated Fourier method.

Example 2.2**.**

Upon the presence of the Dirichlet eigen-elements of the Laplacian operator, we choose

[TABLE]

Recall the Gevrey type class of functions of order $q>0$ and index $p>0$ defined by the spectrum of the Laplacian as follows:

[TABLE]

This is a Hilbert space equipped with the following inner product and norm:

[TABLE]

With the help of the Parseval identity, it is then easy to see that $\bar{\mathbb{W}}=\mathbb{G}_{\overline{M}T,2}$ and $\bar{C}_{0}=\overline{M}$ in Definition 1.3. Now taking $\mathbf{P}_{\varepsilon}^{\beta}:=\overline{M}\Delta+\mathbf{Q}_{\varepsilon}^{\beta}$ to absorb high frequencies in the Laplacian, we get

[TABLE]

*and therefore, it holds that $C_{1}=1/T$ in Definition 1.4. *

Remark 2.3**.**

The conventional cut-off function for the locally Lipschitz $F$ can be taken by

[TABLE]

where $\ell^{\varepsilon}=\ell(\varepsilon)>0$ satisfying $\lim_{\varepsilon\to 0^{+}}\ell^{\varepsilon}=\infty$ is called the cut-off parameter. In the same manner, we have the cut-off function $G_{\ell^{\varepsilon}}(x,t;u;v)$ and thus one can prove that

[TABLE]

*and one can observe that $\lim_{\varepsilon\to 0^{+}}L(\ell^{\varepsilon})=\infty$. The proof can be detailed in [26, 28]. *

3 Main results

3.1 Approximation by the QR method

This part is devoted to establishing an approximate problem of (1) by using the QR method. We begin by adding the perturbing operator to the original system as follows:

[TABLE]

Hereby, we define new diffusion-like coefficients $\overline{\mathcal{D}}_{i}$ for $i=1,2$ in such a way that $\overline{\mathcal{D}}_{i}:=\overline{M}-\mathcal{D}_{i}$. This way one can prove that $\overline{\mathcal{D}}_{i}\in(\overline{M}-M_{1},\overline{M})$ for $M_{1}<\overline{M}$ by $\left(\text{A}_{1}\right)$. Thus, we rely on the choice of the stabilized operator (i.e. $\mathbf{P}_{\varepsilon}^{\beta}:=\overline{M}\Delta+\mathbf{Q}_{\varepsilon}^{\beta}$) to arrive at

[TABLE]

where $F_{\ell^{\varepsilon}}$ and $G_{\ell^{\varepsilon}}$ are defined as in Remark 2.3. This is the coupled regularized system of (1) that we wish to scrutinize in this work. Due to Lemma 2.1, we associate (6) with the Dirichlet boundary condition and the following terminal conditions:

[TABLE]

Hence, these equations form our regularized problem $\left(P^{\varepsilon,n}\right)$.

To study the weak solvability of the regularized problem $\left(P^{\varepsilon,n}\right)$, we use an exponential weight function $e^{\rho_{\varepsilon}(t-T)}$, where $\rho_{\varepsilon}>0$ is called as an $\varepsilon$-dependent auxiliary parameter, to consider the mappings $U^{\varepsilon,n}=e^{\rho_{\varepsilon}(t-T)}u^{\varepsilon,n}$ and $V^{\varepsilon,n}=e^{\rho_{\varepsilon}(t-T)}v^{\varepsilon,n}$. Thus, (6) becomes

[TABLE]

endowed with the Dirichlet boundary condition and with the same terminal data as (7). Henceforward, we define a weak formulation of this transformed system in the following type.

Definition 3.1**.**

For each $\varepsilon>0$, a pair of functions $\left(U^{\varepsilon,n},V^{\varepsilon,n}\right)$ is said to be a weak solution of (8) if

[TABLE]

and it holds that

[TABLE]

*for all $\psi_{1},\psi_{2}\in H_{0}^{1}(\Omega)$. *

It is worth mentioning that to be successful with the Galerkin-type aid we choose in [26], the weight $\rho_{\varepsilon}$ must be large, controlled by the largeness of the magnitude stability of the regularized problem (cf. [24]) and of the cut-off parameter $\ell^{\varepsilon}$. In the following, we only state the well-posedness result of the “scaled” problem (8), while details of proof can be deduced as in [26].

Theorem 3.2**.**

*Assume $\left(\text{A}_{1}\right)$–$\left(\text{A}_{3}\right)$ hold. For each $\varepsilon>0$, the regularized problem (8) admits a pair of weak solutions $\left(U^{\varepsilon,n},V^{\varepsilon,n}\right)$ in the sense of Definition 3.1. Moreover, one has $U^{\varepsilon,n},V^{\varepsilon,n}\in C([0,T];L^{2}(\Omega))$ and $U^{\varepsilon,n}_{t},V^{\varepsilon,n}_{t}\in L^{2}(0,T;(H^{1}(\Omega))^{\prime})$. *

In general, one can get the strong solution of (8) in the sense that $U^{\varepsilon,n}_{t},V^{\varepsilon,n}_{t}$ and $\Delta U^{\varepsilon,n},\Delta V^{\varepsilon,n}$ belong to $L^{2}(0,T;L^{2}(\Omega))$ by increasing the regularity of the corresponding terminal conditions $U^{\varepsilon,n}_{f},V^{\varepsilon,n}_{f}$ (cf. (7)) in $H^{1}(\Omega)$; see [3, Theorem 3.3] for detailed techniques that we can adapt. However, this augment is not practical as in the context we are dealing with the noisy data. In other words, their gradients cannot be measured and even if it is possible, it eventually requires very much effort and expense. Accordingly, this explains why at present, we cannot fully adapt the error estimates for the finite element solution of (8) obtained in, e.g., [9, 6] in Section 4. This open question will be explored in the near future.

3.2 Error analysis

In the following, $n$ will be dependent of $\varepsilon$ (i.e. $n:=n(\varepsilon)\in\mathbb{N}$) due to the argument obtained in Lemma 2.1, showing that $n$ cannot be arbitrarily large. Unlike the previous section where we design the approximate problem, here we denote by $u_{\beta}^{\varepsilon,n},v_{\beta}^{\varepsilon,n}$ the regularized solutions due to the involved regularization parameter $\beta$, recalled from Section 1.3. Note here that the weight $\rho_{\beta}$ we use in proofs of the main results plays the same role as the weight $\rho_{\varepsilon}$ considered in Section 3.1, i.e. it must be large as driven by the smallness of $\varepsilon$. Additionally, we below assume that $C_{1}T\leq 1$ and

[TABLE]

Theorem 3.3** (error estimate for $0<t<T$).**

Assume $\left(\text{A}_{1}\right)$–$\left(\text{A}_{3}\right)$ and $\left(\text{A}_{1}^{\prime}\right)$ hold. Consider the coupled Dirichlet system (1) with terminal data (2) satisfying (3)–(4). Suppose that it has a unique pair of solutions satisfying the source condition

[TABLE]

where the function space $\bar{\mathbb{W}}$ is obtained from the choice of the perturbing operator $\mathbf{Q}_{\varepsilon}^{\beta}$ in Definition 1.3. By the resulting stabilized operator $\mathbf{P}_{\varepsilon}^{\beta}$ in Definition 1.4, we consider $\left(u^{\varepsilon,n}_{\beta},v^{\varepsilon,n}_{\beta}\right)$ as a pair of solutions of the approximate system (6)–(7). Then for $0<\kappa\leq 2C_{1}t$ the following estimate holds:

[TABLE]

*where $C_{3}>0$, $\tilde{C}_{p}>0$ are not dependent of $\varepsilon$. *

There is no doubt that it is hard to attain the convergence of regularization at $t=0$. Herewith, by requiring more information of the source condition we find $t^{\varepsilon}>0$ in such a way that $u_{\beta}^{\varepsilon}(\cdot,t=t^{\varepsilon})$ will be a good approximation of $u(\cdot,t=0)$, involving the expectation operator.

Theorem 3.4** (error estimate for $t=0$).**

Under the assumptions of Theorem 3.3, we further assume that $u,v\in C^{1}(0,T;L^{2}(\Omega))$. Then there exists $t^{\varepsilon}\in(0,T)$ approaching 0 such that

[TABLE]

*where $C_{4}>0$ is independent of $\varepsilon$ and $0<\kappa^{\varepsilon}<\min\left\{2C_{1}t^{\varepsilon},1\right\}$. *

3.2.1 Proof of Theorem 3.3

Put $\mathcal{X}_{\beta}^{\varepsilon,n}\left(x,t\right)=e^{\rho_{\beta}\left(t-T\right)}\left[u_{\beta}^{\varepsilon,n}\left(x,t\right)-u\left(x,t\right)\right]$, $\mathcal{Y}_{\beta}^{\varepsilon,n}\left(x,t\right)=e^{\rho_{\beta}\left(t-T\right)}\left[v_{\beta}^{\varepsilon,n}\left(x,t\right)-v\left(x,t\right)\right]$ for $\rho_{\beta}>0$. Then the “scaled” difference equation for $u$ is given by

[TABLE]

In parallel, we also derive that the difference equation is endowed with the zero Dirichlet boundary condition $\mathcal{X}_{\beta}^{\varepsilon,n}=0$ on the boundary $\partial\Omega$ and with the terminal condition $\mathcal{X}_{\beta}^{\varepsilon,n}\left(x,T\right)=U_{f}^{\varepsilon,n}\left(x\right)-u_{f}\left(x\right)$ for $x\in\Omega$. Hereby, we multiply this equation by $\mathcal{X}_{\beta}^{\varepsilon,n}$ and then integrate the resulting equation over the domain of interest $\Omega$ to arrive at

[TABLE]

Cf. Definition 1.4, we bound $I_{1}$ from below by

[TABLE]

aided by the Cauchy–Schwarz inequality. Meanwhile, we use the Young inequality in combination with Definition 1.3 to get

[TABLE]

Now, by decreasing $\varepsilon$ we can choose the cut-off parameter $\ell^{\varepsilon}$ large in such a way that $\ell^{\varepsilon}\geq\max\left\{\left\|u\right\|_{L^{\infty}\left(0,T;L^{2}\left(\Omega\right)\right)};\left\|v\right\|_{L^{\infty}\left(0,T;L^{2}\left(\Omega\right)\right)}\right\}$. Thus, it holds that $F_{\ell^{\varepsilon}}\left(u;v\right)=F\left(u;v\right)$, according to Remark 2.3. This way we bound $I_{3}$ from below by

[TABLE]

Using $\left(\text{A}_{1}^{\prime}\right)$, the estimate of $I_{4}$ is thus given by

[TABLE]

Grouping the estimates of $I_{i}$ for $i=\overline{1,4}$, we obtain the fact that

[TABLE]

By doing the same token with the difference equation for $v$, we can find that

[TABLE]

Combining (12), (13), taking the expectation and recalling Lemma 2.1, we choose

[TABLE]

and thus, it yields

[TABLE]

with $\tilde{C}_{p}>0$ given by

[TABLE]

By the back-substitutions of $\mathcal{X}_{\beta}^{\varepsilon,n}$ and $\mathcal{Y}_{\beta}^{\varepsilon,n}$, one has

[TABLE]

where we have denoted by $C_{2}\left(\ell^{\varepsilon}\right)=\frac{1}{4}+4L(\ell^{\varepsilon})+\frac{\tilde{L}^{2}\tilde{C}_{p}}{\overline{M}-M_{1}}$

Observe in (14) that if we take $n$ and $\gamma$ satisfying (11) and choose $C_{1}$ such that $C_{1}T\leq 1$, we can go on the estimate (14) as follows:

[TABLE]

It now remains to deal with a fine control of $C_{2}(\ell^{\varepsilon})$ in (15) that involves the blow-up profile of $L(\ell^{\varepsilon})$ (see again in Remark 2.3) as $\varepsilon\to 0^{+}$. In fact, due to the conventional logarithmic rate of convergence in regularization of parabolic problems, our attempt is to bound this $\varepsilon$-dependent constant in such a way that its growth does not destroy the aimed speed. To do that, we basically follow the strategy in [26, Section 5], which enables us to choose

[TABLE]

where $\kappa>0$ is an $\varepsilon$-independent constant. Thus, it gives $e^{8(T-t)L(\ell^{\varepsilon})}\leq\log^{\kappa}(\gamma(T,\beta))$ and consequently, it holds that

[TABLE]

where $C_{3}=\frac{1}{4}+\frac{\tilde{L}^{2}\tilde{C}_{p}}{\overline{M}-M_{1}}>0$ is now independent of $\varepsilon$.

Hence, by choosing $\kappa:=\kappa(t)\leq 2C_{1}t$ we complete the proof of the theorem.

3.2.2 Proof of Theorem 3.4

Using the triangle inequality and aided by Theorem 3.3, one has

[TABLE]

Here, $\kappa^{\varepsilon}:=\kappa(t^{\varepsilon})\in(0,2C_{1}t^{\varepsilon}]$. Observe that since with respect to $t^{\varepsilon}$, the term $\gamma^{-C_{1}t^{\varepsilon}}$ is decreasing, while $t^{\varepsilon}$ itself increases linearly. Thus, for each $\varepsilon>0$ we can find a unique $t^{\varepsilon}\in(0,T)$ such that $t^{\varepsilon}=\gamma^{-C_{1}t^{\varepsilon}}(T,\beta)$ or equivalently, $\frac{\log(t^{\varepsilon})}{t^{\varepsilon}}=-C_{1}\log(\gamma(T,\beta))$. Thanks to the inequality $\log(a)>-a^{-1}$ for all $a>0$, we find that

[TABLE]

By the growth of $\gamma(T,\beta)$ as $\varepsilon\to 0^{+}$, one knows that $\lim_{\varepsilon\to 0^{+}}t^{\varepsilon}=0$. Hence, we complete the proof of the theorem.

3.2.3 Essential remarks

Upon our main results reported in Theorems 3.3 and 3.4, we single out some important remarks where we believe that it is essential to highlight our novel QR method.

• •

By enjoying the exponential weight function in proofs of these theorems, it is clear that the speeds of convergence are exponentially controlled with respect to time and this control becomes larger when our reconstruction runs close to $t=0$.

• •

When $F$ and $G$ are globally Lipschitz-continuous with respect to the arguments $u$ and $v$, we obtain the Hölder-type rate for $t\in(0,T)$, which means $\mathcal{O}\left(\gamma^{-2C_{1}t}(T,\beta)\right)$, whilst it is of the logarithmic order $\mathcal{O}(\log^{-1}(\gamma(T,\beta)))$ at $t=0$. These are what we have obtained in [26, Theorem 4.5]. The presence of the term $\log^{\kappa}(\gamma(T,\beta))$ in Theorem 3.3 shows that the locally Lipschitz-continuous nonlinearities usually encountered in real-world applications require much effort to approximate. In the near future, we will figure out a better way to design a new bound for the cut-off constant $L(\ell^{\varepsilon})$ (cf. (16)).

• •

Cf. [26], if the model under consideration involves the gradient terms accounting for, e.g., haptotaxis in a population of cells (as presented in a model for haptotaxis and chemotaxis effects on cells’ motion; cf. [27]), our proposed method is still valid. A slightly lower rate of convergence is expected in controlling the gradient terms, but theoretically, it still converges logarithmically.

• •

If the condition (11) is changed to

[TABLE]

we eventually obtain the convergence rate $\mathcal{O}\left(\varepsilon^{1/2}+\log^{-1}(\gamma(T,\beta))\right)$ in the space $C(t,T;L^{2}(\Omega))\cap L^{2}(t,T;H^{1}_{0}(\Omega))$ for any $0\leq t<T$. This is easily deduced from the estimate (14). Therefore, the presence of Theorem 3.4 as well as the additional regularity assumption that $u,v\in C^{1}(0,T;L^{2}(\Omega))$ can be neglected. This is also a new finding compared to the convergence result in our original work [26].

• •

The special difference between the deterministic noise in [26, 24, 18, 28, 22] and the stochastic noise we are taking into account is the following. It is not fully about the involved expectation operator, but about the findings of $K_{0}$ and $K_{1}$ in (11). It shows that our QR method in this context really needs a very careful adaptation because not only $n$ has to be chosen properly, but also the eigenvalues driven by such $n$ need a fine control with respect to $\varepsilon$ and $\gamma$. Consider $\Omega$, for example, as an open parallelepiped $(0,a_{1})\times\ldots(0,a_{d})\subset\mathbb{R}^{d}$ with $a_{i}>0$, $i\in\left\{1,\ldots,d\right\}$. As one of advantages of this QR method, we only need to solve the simple Dirichlet eigenvalue problem regardless of the complex structure involved in the diffusion term, whenever it is essentially bounded. Henceforth, the Dirichlet eigen-elements are given by

[TABLE]

In this circumstance, we mean $n$ in (11) as $\left|n\right|$ where $n:=\left(n_{1},\ldots,n_{d}\right)$. Then, we can choose $n_{j}\propto\varepsilon^{-2\theta}$ for $j=1,\ldots,d$ and $\theta\in(0,1)$. Notice that $\mu_{n}\propto\left|n\right|^{2}$, then $\gamma^{C_{1}T}(T,\beta)\propto\varepsilon^{\max\left\{\theta-1,-4\theta p\right\}}$ is needed. Essentially, we have drawn a possibility to choose $n$ and $\gamma$ in the context of (11). Additionally, one can take $\beta=\varepsilon$ and thus $\gamma(t,\beta)=\varepsilon^{\frac{t}{C_{1}T^{2}}\max\left\{\theta-1,-4\theta p\right\}}$. Hence, the following error bounds hold

[TABLE]

[TABLE]

where $\kappa(t)\in(0,2C_{1}t]$ and $\kappa^{\varepsilon}\in\left(0,\min\left\{2C_{1}t^{\varepsilon},1\right\}\right)$.

• •

In case of deterministic noises, the current analysis only needs the true terminal data in $L^{2}(\Omega)$. Moreover, the assumptions (11) will be as we required in [26, Theorem 4.5].

4 Numerical tests

Throughout the examples, we choose $\Omega=(a,b)=(0,\pi)$ as a sample domain of tissue in one-dimensional. Accordingly, the Dirichlet eigen-elements are elementary and have the following form:

[TABLE]

Given $\varepsilon>0$, as discussed in Section 3.2.3 we choose that $n=\left\lfloor\varepsilon^{-2\theta}\right\rfloor$ for $\theta\in(0,1)$, $\beta(\varepsilon)=\varepsilon$ and $\gamma(t,\varepsilon)=\varepsilon^{\frac{t}{C_{1}T^{2}}\max\left\{\theta-1,-4\theta p\right\}}$ for all $t\in[0,T]$. Then (11) is satisfied with $K_{0},K_{1}\leq 1$. Employing Example 2.2, we define the stabilized operator $\mathbf{P}_{\varepsilon}^{\beta}$ as follows:

[TABLE]

where $\mathbb{B}_{\varepsilon}:=\left\{j\in\mathbb{N}:\mu_{j}\leq\log\left(\varepsilon^{\frac{1}{C_{1}T}\max\left\{\theta-1,-4\theta p\right\}}\right)\right\}$ is the admissible set of our Fourier frequencies, and $C_{1}=\overline{M}$.

We generate the discrete observations $u_{f}^{\varepsilon,j},v_{f}^{\varepsilon,j}$ for $1\leq j\leq n$ by the assumptions (3)–(4) where $\left\langle\xi_{i},\phi_{j}\right\rangle$ for $i=1,2$ are random variables with normal distribution. In our illustration, we employ the manufactured solution to verify the validity of the regularization scheme. When doing so, we know the explicit expression of the true $u_{f}$ and $v_{f}$. Cf. [16, Section 4], the adaptive Filon-type method is rather reliable to compute the inner products $\left\langle u_{f},\phi_{j}\right\rangle$ and $\left\langle v_{f},\phi_{j}\right\rangle$ due to the low level of frequency restricted by the stabilized operator. It then enables us to obtain the reconstructions $U_{f}^{\varepsilon,j},V_{f}^{\varepsilon,j}$.

We now take into account the numerical regularized solution by using the backward Euler method. In this regard, a uniform grid of mesh-points $(x,t)=(x_{m},t_{k})$ is used. Here $x_{m}=a+m\Delta x$ and $t_{k}=k\Delta t$ for $0\leq m\leq M$ and $0\leq k\leq K$ where $\Delta x=\frac{b-a}{M}$, $\Delta t=\frac{T}{K}$ and $M,K\in\mathbb{N}^{*}$. Henceforth, we find $u_{m,k}^{\varepsilon,n}\approx u^{\varepsilon,n}(x=x_{m},t=t_{k})$ and $v_{m,k}^{\varepsilon,n}\approx v^{\varepsilon,n}(x=x_{m},t=t_{k})$ by solving the following fully discretized system:

[TABLE]

for $1\leq m\leq M-1$ and $0\leq k\leq K-1$. Accordingly, we endow this system with the Dirichlet conditions $u^{\varepsilon,n}_{0,k}=u^{\varepsilon,n}_{M,k}=v^{\varepsilon,n}_{0,k}=v^{\varepsilon,n}_{M,k}=0$ and the terminal conditions $u^{\varepsilon,n}_{m,K}=U^{\varepsilon,n}_{f}(x=x_{m})$, $v^{\varepsilon,n}_{m,K}=V^{\varepsilon,n}_{f}(x=x_{m})$. Note that our objective here is to solve this system backwards in time and thus, $k$ runs from $K$ to [math] in the iterations. Then, for each $k=0,\ldots,K-1$ we arrive at a linear algebraic system with a Toeplitz-like matrix, as follows:

[TABLE]

where the solution $\mathbb{W}^{\varepsilon,n}_{k}\in\mathbb{R}^{2(M-1)}$, the block matrix $\mathbb{A}\in\mathbb{R}^{(2M-2)\times(2M-2)}$ and the function $\mathbb{F}_{\varepsilon}\in\mathbb{R}^{2(M-1)}$ are defined by

[TABLE]

[TABLE]

[TABLE]

with $\mathbb{D}_{1}^{k+1},\mathbb{D}_{2}^{k+1}\in\mathbb{R}^{K-1}$ taken by

[TABLE]

[TABLE]

for $\overline{\mathcal{D}}_{1}^{k+1}=\overline{\mathcal{D}}_{1}\left(u_{k+1}^{\varepsilon,n}\right)$, $\overline{\mathcal{D}}_{2}^{k+1}=\overline{\mathcal{D}}_{2}\left(v_{k+1}^{\varepsilon,n}\right)$ and $\bar{\alpha}=\frac{\Delta t}{\Delta x^{2}}$.

We remark that for each iteration step $k$, (18) has a solution because the block matrix $\mathbb{A}$ is symmetric positive definite. In this illustration, we only seek the approximation of $u^{\text{ex}}(x_{m},t=t^{\varepsilon})$, but usually $t^{\varepsilon}$ does not belong to the underlying mesh-points. We then rely on the nearest point of this $t^{\varepsilon}$, which means that for each $\varepsilon>0$ the error, denoted by $E(t_{k})$, is computed for $t_{k}$ such that $t_{k}\geq t^{\varepsilon}$ and $t_{k-1}<t^{\varepsilon}$. Given the number of samples $\mathcal{M}\in\mathbb{N}^{*}$ and choosing $T=1$, the averaged error $E(t_{k})$ reads as

[TABLE]

where $u^{\text{ex}}$ denotes the manufactured solution that we will choose in the numerical tests below. Note that for simplicity, we eventually omit the index $m_{s}$ in the regularized solution in the above numerical setting. The value of $t^{\varepsilon}$ can be approximately computed by the equation $\log(t^{\varepsilon})=-t^{\varepsilon}\log\left(\varepsilon^{\frac{1}{T}\max\left\{\theta-1,-4\theta p\right\}}\right)$ (cf. Section 3.2.2).

4.1 Test 1

We begin the numerical verification by the most widely used Lotka–Volterra model for tumor growth kinetics where in (1) one has $\mathcal{D}_{i}$ are positive constants and the corresponding Lotka–Volterra functions are given by

[TABLE]

Here, the involved parameters in $F,G$ are positive constants and cf. [10], they possess the following biological meaning: $b_{u}$ is the negative effects of tumor on normal cells, $b_{v}$ is a lumped, phenomenological term presented by, e.g., the immune response; $r_{u},r_{v}$ are, respectively, maximum growth rates of normal and tumor cells; $K_{u},K_{v}$ denote their maximal densities.

Given $\mathcal{D}_{1}=0.5$ and $\mathcal{D}_{2}=1$ for slow diffusions, which then lead to $\overline{M}=1.1$ and $M_{1}=1$, we choose $r_{u}=r_{v}=K_{u}=K_{v}=1$ and $b_{u}=2,b_{v}=0.5$, usually describing the slight domination of tumor cells. The manufactured solutions are given by $u^{\text{ex}}=e^{-t}\sin(x)$ and $v^{\text{ex}}=t^{2}x(\pi-x)$, ensuring the positivity of solutions. This way we need to take

[TABLE]

By the choice of $F$ and $G$, we find that $L(\ell^{\varepsilon})=1+3\ell^{\varepsilon}$ and hence, for every time step $\ell^{\varepsilon}$ can be controlled by

[TABLE]

where we choose $\kappa(t)=C_{1}t>0$. Thus, we can work with the cut-off functions $F_{\ell^{\varepsilon}}$ and $G_{\ell^{\varepsilon}}$ in this circumstance.

4.2 Test 2

In the second example, we take into account the nonlocal problem as in [7], where $F,G$ are independent of the solutions, i.e. $F(x,t;u;v)=F(x,t)$ and $G(x,t;u;v)=G(x,t)$. Then the diffusion terms are particularly expressed by

[TABLE]

Basically, this context can be understood that the tumor cells are on the way to dominate the normal cells, since the evolution system is now uncoupled and partly driven by the total population of each species. In this test, we choose $u^{\text{ex}}=0.25\log(3+t)\sin^{2}(x)$ and $v^{\text{ex}}=t\sin(x)$ with $d_{1}=0.01$ and $d_{2}=0.05$. Then, it follows that $\overline{M}=4,M_{1}=3.5$ and

[TABLE]

Note that due to the nonlocal diffusion, we compute $\overline{\mathcal{D}}_{i}^{k+1}$ ($i=1,2$) involved in the algebraic system (18) by employing the midpoint rule. It yields

[TABLE]

where $\tau\left(u\right)=\left|u\left(1-u\right)\right|$. It is worth noting that the choice of the numerical integration for the diffusion term has to be careful. Although the efficient Gauss–Legendre quadrature and its variants seem applicable by choosing suitable nodes, in practice one cannot get measurements exactly at those abscissae. One can also take into account the Newton–Cotes method for the uniform grid in this implementation, but this approach is usually unstable when $M$ is large. Therefore, one can see the impediment in handling regularization methods in numerics, especially when there are numerous nonlinear factors involved in the model.

4.3 Numerical results

We take $p=1$ and $\theta=0.3$ in these two tests. The choice of $\theta$ indicates that we aim to see the convergence for a not very large number of steps of discrete observations $n$. Meanwhile, the value of $\varepsilon$ varies from $10^{-3}$ to $10^{-5}$.

Our convergence results are concluded in Figs. 1 and 2, where we have computed the numerical errors of two tests with $\mathcal{M}=100$ samples (the number of samples is implemented in the horizontal line of each figure) and $M=15$, $K=100$. Having these errors allows us to obtain the averaged errors corresponding to each solution, which afterward confirms our theoretical expectation. Particularly, in Test 1 the averaged error for $u$ decreases from $0.01224$ (for $\varepsilon=10^{-3}$) to $0.00284$ (for $\varepsilon=10^{-5}$), while for $v$ it slowly goes down to $0.03971$ from $0.06545$. Similarly, in Test 2 the averaged error for $u$ is from $0.00191$ (for $\varepsilon=10^{-3}$) to $0.00177$ (for $\varepsilon=10^{-5}$), while for $v$ it decreases from $0.04640$ to $0.02689$. Henceforth, one can see that our approximation is acceptable from $\varepsilon=10^{-3}$, i.e. this time the regularized solution is close to the true solution.

5 Conclusions

We have improved and adapted the recently developed QR method in [26] to regularization of a terminal-boundary value parabolic problem with white Gaussian noise. Although this noise process is very standard in the stochastic setting, this attempt concretely consolidates the mathematical quality of the reconstruction method we have established so far, where both the well-posedness of the regularized problem and the convergence analysis are successfully obtained. Besides, it again confirms the conventional logarithmic rate that we usually meet in regularization of parabolic problems. This time we introduce a new type of perturbation that allows us to compute the stabilized operator in the numerical verification of the scheme.

Although several questions remain open in [26], we would like to add one more interesting problem. In view of the localization of the brain tumor source in this paper, one can address the model in a domain with moving boundaries as investigated in, e.g., [2]. This way may be helpful in brain pathologies, especially in improving the Magnetic Resonance Imaging method as one wants to anatomically see the motion of the brain (see, e.g., [12]).

Acknowledgments

V.A.K thanks Thi Kim Thoa Thieu (Gran Sasso Science Institute, Italy and Karlstad University, Sweden) for fruitful discussions during the completion of the manuscript. V.A.K thanks Prof. Loc Nguyen (North Carolina, USA) for the hospitality he provided.