Numerical differentiation for two-dimensional scattered data on arbitrary domain base on Hermite extension with an implicit iteration process

1.   Introduction

In the application of many mathematical physics problems, we need to estimate derivatives of an unknown function from given noisy data. It turns out to be an ill-posed problem, which means, the small errors in the measurement data can induce huge errors in its computed derivatives. Many methods and techniques have been proposed regarding this topic ^{[1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17]}. According to type of regularization techniques, these methods can be classified into finite difference methods, mollification methods, differentiation by integration method and Tikhonov methods. Most of these methods are for one-dimensional case, but only a few are for high-dimensional cases ^[18,19,20]. These methods of dealing with two-dimensional problems are basically aimed at the regular region, and there are few methods that can deal with the numerical differentiation problem on the irregular domain. For example, if the data are given at scatter points of the domain in Figure 1, numerical implementation of most existing methods is difficult.

By the extension theorem ^[21], any function on sub-domain of $\mathbb{R}^2$ can be extended to $\mathbb{R}^2$ without losing its smoothness. There are many work of extension problem has been developed since Whitney's seminal work ^{[22,23,24,25]}. And it is well known that functions with high smoothness can be approximated very precisely by their Hermite expansion. So in this paper, we consider a Hermite extension method to deal with the numerical differentiation problem on the irregular domain. That is to say, we take a function defined on any domain as part of a function on $\mathbb{R}^2$. For the one-dimensional case, the numerical differentiation method along the line of this method has been given in ^[26]. The numerical process of extension is usually unstable, so regularization technology is needed. As an alternative regularization method, we use a modified implicit iteration method in this paper. Compared with Tikhonov method which is used in ^[17], the implicit iteration method can select larger regularization parameters in each iteration, which makes the calculation process more stable. In ^[27], Jin has given the detailed theory and numerical implementation of implicit iteration method in Hilbert scales. But the application of the method to our problem still needs some improvement. The existing implicit iterative method can only deal with the problem of finite smoothness, and the numerical implementation of the algorithm is inconsistent for different smoothness. In this paper, we will present a modified form of implicit iteration method. It can deal with numerical differentiation of functions with any smoothness and the solution process is uniform.

This paper is organized as follows: In the next section, we introduce some preliminary materials. In Section 3, we describe the modified implicit iteration method for numerical differentiation problem and give some auxiliary results. The convergence estimate of the approximation solution can be founded in Section 4. In Section 5, some numerical examples are given. Conclusion is then given in the final Section 6.

2.   Preliminaries

In this section, we introduce some notations and preliminaries that will be used throughout the paper.

2.1.   Hermite functions in $\mathbb{R}^2$

Let ${\bf{x}} = (x_1, x_2)$, $|{\bf{x}}| = \sqrt{x_1^2+x_2^2}$. Let

be the 2-dimensional Fourier transform of the function $f({{\bf{x}}})\in L^2(\mathbb{R}^2)$. The corresponding inverse Fourier transform of the function $\hat{f}(\mathit{\boldsymbol{\omega}})$ is

And $\|\cdot\|_p$ denotes the norm of the Sobolev space $H^p(\mathbb{R}^2)$ with $p\geq 0$ defined by

Particularly, for $p = 0$ we can recover the $L^2(\mathbb{R}^2)$ norm, i.e.,

Let ${\mathit{\boldsymbol{\ell}} = (\ell_1, \ell_2)}$, ${\ell}_i (i = 1, 2)$ being non-negative integers. Set ${|{{\mathit{\boldsymbol{ \ell }}}}|_1 = \ell_1+\ell_2}$ and ${\mathit{\boldsymbol{\ell}}\cdot{\bf{x}} = \ell_1x_1+\ell_2x_2}$. The normalized 2-dimensional Hermite function is defined by

where

We know that the Fourier transform of the Hermite functions can be given as^[28]:

The set of Hermite functions satisfy the orthogonality relations

The Hermite expansion of a function $f\in L^2(\mathbb{R}^2)$ is as

with the Fourier-Hermite coefficients

2.2.   Problem description

Suppose that $\Lambda$ is a subdomain of $\mathbb{R}^2$. For any two-tuples ${\alpha} = (\alpha_1, \alpha_2)\in \mathbb{N}^2, |\alpha|_1 = \alpha_1+\alpha_2$. The notation $\partial_j$ stands for $\frac{\partial}{\partial x_j}$ and $D^{\alpha}f = \partial^{\alpha_1}_{x_1}\partial^{\alpha_2}_{x_2}f$. For any positive integer $q$, we define $D^qf = :\{D^{\alpha}f:|\alpha|_1 = q\}$ and

The norm $\|\cdot\|_{s, \Lambda}$ in Sobolev space $H^s(\Lambda)$ is defined as

where $s = 0$ and $\|\cdot\|_{0, \Lambda}$ denotes the $L^2(\Lambda)$ norm.

Suppose that $g(x)\in H^{p}(\mathbb{R}), p\geq 2$ and we only know its approximate function $g^{\delta}$ on $\Lambda$ such that

where $\delta > 0$ is a given constant called the error level. Our problem is to calculate approximate derivatives of $g$ on $\Lambda$ from the noisy data $g^{\delta}$, or, equivalently, to construct a function $f^{\delta}(x)$ from $g^{\delta}(x)$ which is close to $g(x)$ in the sense that

3.   Implicit iteration method for numerical differentiation

For any vector $\vec{\mathsf{f}} = \{\mathsf{f}_{\mathit{\boldsymbol{\ell}}}\}_{|{{\ell}}|_1 = 0}^{\infty}\in l^2$, if we let

then the process of constructing an approximation function $f^{\delta}$ from data $g^{\delta}$ can be transformed to solving the following equations

In this paper, we present an modified implicit iteration method to solve the above equations. For this purpose, we introduce the following operator:

It is obvious that $\mathcal{R}$ is unbounded self-adjoint strictly positive definite operator. Then we choose

as the approximation of $g$, where $\vec{\mathsf{f}}^{\delta}_n$ is determined by the following implicit iteration process

where $\beta_k > 0$ are properly chosen real numbers. For reference ^[27], the positive number

plays the role of the regularization parameter and we will chosen it as the solution of the nonlinear equation

with a constant $C\geq 1$. If we let $\mathcal{T} = \mathcal{H}\mathcal{R}^{-1}$, then $\vec{\mathsf{f}}^{\delta}_{n}$ possesses the representation ^[27]

Remark 3.1. If we use the operator $\mathcal{B} = \left(\sum_{|\alpha|_1 = q}D^{\alpha}\right)$ with some a constant $q$ instead of $\mathcal{R}$, then we return to the framework in ^[27] and the convergence results can be obtained accordingly. When $p$ and $q$ satisfy a certain relation, the result is order optimal. It should be noticed that for large $q$, the numerical process of the method is difficult. We will point that the method is always order optimal when we use the operator $\mathcal{R}$ and the numerical process is uniform for any $p$.

The following lemma holds for $s_n(\lambda)$.

Lemma 3.1. ^[27] The function $s_n:(0, c]\rightarrow (0, \infty)$ with $c = \|\mathcal{T}\|^2$ and the corresponding residual function $r_n(\lambda): = 1-\lambda s_n(\lambda)$ obey the properties

From above lemma, we can deduce the following results.

Lemma 3.2.

Proof. For $\lambda\leq \sigma^{-1}$,

and

Moreover, for $\lambda\geq\sigma^{-1}$,

and

4.   Error estimate

Owing to $g\in H^p(\mathbb{R})$, we suppose that

where $E$ is a constant. Set the vector $\vec{\mathsf{g}}$ contains all Fourier-Hermite coefficients of $g$, i.e.,

Let

We define the vector $\vec{\mathsf{f}}_n$ as

then we have

In our further analysis, we shall make use of the following lemmas.

Lemma 4.1. If the condition (4.1) holds, then

where

Proof. From (2.1) and (3.3), we can obtain

and

Lemma 4.2. If the condition (4.1) holds, we have

and

Proof. From (2.10), (3.7), (4.10) and the triangle inequality

And by using the triangle inequality, (2.10), (3.10) and (4.8)–(4.10)

Moreover, in terms of the triangle inequality, (2.10), (3.9), (3.10) and (4.7), we have

Lemma 4.3. ^[21] Let $\Omega$ be a domain in $\mathbb{R}^2$ satisfying the cone condition. There exists a constant $K$ depending on $\epsilon_0$ and $j$, $s$, such that for any $0 < \epsilon\leq \epsilon_0$ and $0\leq j\leq s$

Lemma 4.4. Suppose that the vector sequence $\vec{\mathsf{h}}_n^{\delta} = \left\{\mathsf{h}^{\delta}_{\mathit{\boldsymbol{\ell}}}\right\}_{|{{\mathit{\boldsymbol{ \ell }}}|_1 = 0}}^{\infty}$ satisfies

then for any $\Omega\subseteq\Lambda$ satisfying the cone condition, there exists a constant $M$

Proof. It is easy to deduce that there exist a constant $\delta_0$ such that

And for simplicity, we prove the theorem with $\delta < \delta_0$. Let

and we have

By Parseval's formula, we can see that the second term $I_2$ satisfies

Hence

So all we need is to prove there exist a constant $M_1$ such that

Note that

and

Therefore

Now we prove (4.28) by using reduction to absurdity, if (4.28) does not hold, then for any $q > p$ there exist a sequence $\delta_i$ such that

If not, $\exists \bar{q}$ for any $\delta\rightarrow 0$

then (4.28) can be derived by 4.29 and Lemma 4.3 with $\epsilon = \left(\frac{\delta_i}{k_3^{p}}\right)^{\frac{\bar{q}-p}{p}}$, $s = \bar{q}$ and $j = p$. Then

Therefore

So there exists a $\bar{\delta}$ such that

which contradicts the assumptions of the Lemma.

Theorem 4.1. Suppose that the conditions (2.10) and (4.1) hold, $f^{\delta}_n$ is defined by (3.4) and (3.7) then for any $\Omega\subseteq\Lambda$ satisfying the cone condition and $0 < j\leq p$,

Proof. Let

Then from Lemma 4.2, we have

Thus, by using Lemma 4.4, there exists a constant $M$

Then

Moreover, by using (2.10), (3.7) and the triangle inequality

The assertion of theorem follows from (4.42), (4.43) and Lemma 4.3.

5.   Numerical realization

The data are usually given at scatter points in practical applications. Let ${\bf{x}}_i\in \Lambda (i = 1, 2, \ldots, m)$ be the given points and

be the noisy data vector. Let ${{\mathit{\boldsymbol{ \sigma }}}}^{({\mathit{\boldsymbol{j}}})} = (\sigma^{(j_1)}, \sigma^{(j_2)}), (0\leq j_1, j_2\leq n)$ being the Hermite-Gauss type interpolation points and ${\mathit{\boldsymbol{ \rho }}}^{({{\mathit{\boldsymbol{j}}}})}$ are the corresponding Hermite-Gauss weights. For $f, g\in L^2(\mathbb{R}^2)$, we define the discrete inner product

and

Let

then we give the matrices $\mathsf{A}_{(n+1)^2\times(n+1)^2}$, $\mathsf{R}_{(n+1)^2\times(n+1)^2}$, $\mathsf{H}_{m\times(n+1)^2}$ as

With these preparations, the discrete form of the implicit iteration method can be given as

Suppose that

Similar to what is done in ^[27], we take $\beta_1 = 1, \beta_k = q^{k-1}\beta_1$ with some $q < 1$ and choose $n$ as the first integer for which

and then adjust the parameter $\beta_n$ such that

where $C_1, C_2$ are two constants that obey $1\leq C_1\leq C_2$.

Remark 5.1. It should be noted that the Hermite-Gauss points are only used to calculate the matrix $\mathsf{R}$, regardless of the location of the noisy data.

6.   Numerical tests

In this section, we give some numerical tests to verify the effect of the new method. All tests are realized on Windows 10 system with Memory 16GB, CPU Intel(R) Core(TM)i7-8550U by using Matlab 2017b. Let $\mathsf{x} = ({\bf{x}}_1, {\bf{x}}_2, \ldots, {\bf{x}}_m)^T$ and the perturbed data are generated by

where $\epsilon$ is the error level and ${\rm{randn}}({\rm{size}}(\cdot))$ is Matlab functions. In all cases we choose the parameter $n = 64$, $q = 1/2$ and $C = 1.01$. We have tested these parameters with other values, and the results are similar. In order to adapt to the characteristic of Hermite function approximation, the scaling factor ^[29] is used in the numerical processing.

Example 6.1. ^[18] Let $\Lambda = \{(x_1, x_2)\; |\; x_1^2+x_2^2\leq 1\}$ is a disk and scatter nodes are given as Figure 2a. We choose the exact function as $g({{\bf{x}}}) = (x_1^2+x_2^2-2)^3$ and set $\epsilon = 0.01$. The numerical results are exhibited in Figures 2c–3h.

Example 6.2. ^[18] Let $\Lambda = \{(x_1, x_2)\; |\; 0\leq x_1\leq 1, \; 0\leq x_2\leq \pi\}$ is a rectangle and scatter nodes are given as Figure 3i. We choose the exact function as $g({{\bf{x}}}) = (x_1^2+x_2^2-2)^3$ and set $\epsilon = 0.01$. We have shown the numerical results in Figures 3k–4h.

Example 6.3. ^[19] Let $\Lambda = \{(x_1, x_2)\; |\; -1\leq x_1\leq 3, \; -1\leq x_2\leq 3\}$, the data are given at the equidistant nodes whose sampling step was $0.1 \times 0.1$. The exact function is chosen as $g({{\bf{x}}}) = \sin\left(\frac{1}{2}x_1^2+\frac{1}{4}x_2^2+3\right)\cos(2x_1+1-\exp(x_2))$. We have shown the numerical results in Figures 4i–5h with $\epsilon = 0.01$.

Example 6.4. Now we let $\Lambda$ is a irregular domain and scatter nodes are given as Figure 5i. We choose the exact function as $g({{\bf{x}}}) = \cos(x_1\cdot x_2)$ and set $\epsilon = 0.01$. The numerical approximations and corresponding errors are shown in Figures 5k–6h.

All the above numerical results show that the proposed method is effective.

7.   Conclusions

In this paper, we present a Hermite extension method with an implicit iteration process for numerical differentiation of two-dimensional functions. Because the method can directly deal with the data given on any domain, it is more convenient than other methods in practical application. The theoretical results show that the convergence rates of the method is self-adaptive.

Conflict of interest

The authors declared that they have no conflicts of interest to this work. We declare that we do not have any commercial or associative interest that represents a conflict of interest in connection with the work submitted.