Galerkin Method with Trigonometric Basis on Stable Numerical Differentiation

TL;DR

This paper develops a Galerkin method with trigonometric basis for stable numerical differentiation of functions on (0, 2π), providing optimal convergence rates, noise-independent parameter choice, and effective handling of discontinuities through numerical examples.

Contribution

The paper introduces a novel Galerkin approach using trigonometric basis for numerical differentiation, achieving optimal convergence and robustness to noise, with explicit parameter selection strategies.

Findings

Achieves optimal convergence rate of O(δ^{2μ/(2μ+1)})

Provides noise-independent parameter choice for specific function forms

Demonstrates effective differentiation of functions with discontinuities

Abstract

This paper considers the $ p $ ($ p=1,2,3 $) order numerical differentiation on function $ y $ in $ (0,2\pi) $. They are transformed into corresponding Fredholm integral equation of the first kind. Computational schemes with analytic solution formulas are designed using Galerkin method on trigonometric basis. Convergence and divergence are all analysed in Corollaries 5.1, 5.2, and a-priori error estimate is uniformly obtained in Theorem 6.1, 7.1, 7.2. Therefore, the algorithm achieves the optimal convergence rate $ O( \delta^{\frac{2\mu}{2\mu+1}} ) \ (\mu = \frac{1}{2} \ \textrm{or} \ 1)$ with periodic Sobolev source condition of order $ 2\mu p $. Besides, we indicate a noise-independent a-priori parameter choice when the function $ y $ possesses the form of \begin{equation*} \sum^{p-1}_{k=0} a_k t^k + \sum^{N_1}_{k=1} b_k \cos k t + \sum^{N_2}_{k=1} c_k \sin k t, \ b_{N_1}, c_{N_2} \neq 0, \end{equation*} In particular, in numerical differentiations for functions above, good filtering effect (error approaches 0) is displayed with corresponding parameter choice. In addition, several numerical examples are given to show that even derivatives with discontinuity can be recovered well.

Tables

Table 1. Table 1: Above experiments correspond to the three examples in Example 8.1. Case I,II are uniformly set as ( δ , δ i ) = ( 0.01 , 0 ) 𝛿 subscript 𝛿 𝑖 0.01 0 (\delta,\delta_{i})=(0.01,0) and ( δ , δ i ) = ( 0.01 , 0.01 ) 𝛿 subscript 𝛿 𝑖 0.01 0.01 (\delta,\delta_{i})=(0.01,0.01) respectively. Notice that r 𝑟 r denotes the relative error.

		$n$	$2$	$4$	$6$	$8$	12
$p=1$	${\delta }_i=0$	$r$	0.4023	0.2132	$1.5194e^{-16}$	$1.5194e^{-16}$	0.0554
	${\delta }_i=0.01$	$r$	0.4024	0.2035	$0.0133$	$0.0152$	0.0532
$p=2$	${\delta }_i=0$	$r$	1.0695	0.6808	$7.5244e^{-15}$	$1.1876e^{-14}$	0.3666
	${\delta }_i=0.01$	$r$	1.0830	0.7046	$0.0732$	0.1051	0.3046
$p=3$	${\delta }_i=0$	$r$	1.1879	1.3802	$6.6497e^{-14}$	$1.6561e^{-13}$	0.0469
	${\delta }_i=0.01$	$r$	1.1884	1.3827	$0.0081$	$0.0150$	0.0355

Table 2. Table 2: t 𝑡 t denotes the CPU time ( s ) 𝑠 (s) for the corresponding experiment in Table 1.

		$n$	$2$	$4$	$6$	$8$	12
$p=1$	${\delta }_i=0$	$t$	0.3341	0.3673	0.4178	0.4577	0.5983
	${\delta }_i=0.01$	$t$	0.3391	0.3728	0.4488	0.4935	0.5748
$p=2$	${\delta }_i=0$	$t$	0.4670	0.6070	0.6385	0.8051	0.8400
	${\delta }_i=0.01$	$t$	0.4892	0.5641	0.6218	0.7126	0.8160
$p=3$	${\delta }_i=0$	$t$	0.2755	0.3958	0.4497	0.5024	0.6810
	${\delta }_i=0.01$	$t$	0.4723	0.5151	0.6348	0.7240	1.1838

Table 3. Table 3: Above experiments correspond to Example 8.2, 8.3 and 8.4 respectively. Case I,II are uniformly set as ( δ , δ i ) = ( 0.01 , 0 ) 𝛿 subscript 𝛿 𝑖 0.01 0 (\delta,\delta_{i})=(0.01,0) and ( δ , δ i ) = ( 0.01 , 0.01 ) 𝛿 subscript 𝛿 𝑖 0.01 0.01 (\delta,\delta_{i})=(0.01,0.01) respectively. Notice that r 𝑟 r denotes the relative error. ( k = 8 𝑘 8 k=8 )

		$n$	$4$	$6$	$8$	$16$	24
$p=1$	${\delta }_i=0$	$r$	0.2786	0.2551	0.2294	0.1474	0.1294
	${\delta }_i=0.01$	$r$	0.2734	0.2486	$0.2216$	$0.1378$	0.1187
$p=2$	${\delta }_i=0$	$r$	0.4148	0.3175	0.2754	0.2068	0.1636
	${\delta }_i=0.01$	$r$	0.4239	0.3323	0.2948	0.2603	0.2667
$p=3$	${\delta }_i=0$	$r$	0.1413	0.1185	0.1209	0.1137	0.1490
	${\delta }_i=0.01$	$r$	0.1446	0.1185	0.1060	0.1383	0.4010

Table 4. Table 4: t 𝑡 t denotes the CPU time ( s ) 𝑠 (s) for the corresponding experiments in Table 3.

		$n$	$4$	$6$	$8$	$16$	24
$p=1$	${\delta }_i=0$	$t$	1.73	2.67	3.89	12.47	103.05
	${\delta }_i=0.01$	$t$	1.55	2.66	4.29	13.80	142.73
$p=2$	${\delta }_i=0$	$t$	3.77	6.57	11.38	73.04	175.93
	${\delta }_i=0.01$	$t$	3.95	6.46	11.09	116.55	230.05
$p=3$	${\delta }_i=0$	$t$	6.10	13.06	19.47	127.01	304.70
	${\delta }_i=0.01$	$t$	7.79	11.45	17.89	149.72	438.42

Equations

$$\sum _{k=0}^{p-1}{a}_k{t}^k+\sum _{k=1}^{N_1}{b}_k\cos kt+\sum _{k=1}^{N_2}{c}_k\sin kt,b_{N_1},c_{N_2}\ne 0,$$

$$\Psi [f]:=\frac{1}{n-1}\sum _{i=1}^{n-1}({\tilde{y}}_i-f(x_i){)}^2+\alpha \Vert f^{\prime \prime }{\Vert }^2$$

$$A^{(p)}\phi :=\frac{1}{(p-1)!}{\int }_0^x(x-t{)}^{p-1}\phi (t)dt=y(x),x\in (0,2\pi ).$$

$$y^{\delta },y,\delta y\in L^2(0,2\pi )\text{and}\Vert \delta y{\Vert }_{L^2}\le \delta ,$$

$$X_n:=span\{\frac{1}{\sqrt{2\pi }},\frac{\cos t}{\sqrt{\pi }},\frac{\sin t}{\sqrt{\pi }},\cdots ,\frac{\cos nt}{\sqrt{\pi }},\frac{\sin nt}{\sqrt{\pi }}\},$$

$$A_n^{(p)}{\phi }_n=y_n{\phi }_n\in X_n,y_n:=P_{n}y\in X_n,$$

$${A_n^{(p)}}^{†}{P}_n{y}^{\delta }\text{with}{n}^{(p)}=n^{(p)}(\delta ),$$

$$n^{(p)}:=n^{(p)}(\delta )\to +\infty ,\delta \to 0^{+}$$

$$\Vert {A_{n^{(p)}(\delta )}^{(p)}}^{†}{P}_{n^{(p)}(\delta )}{y}^{\delta }-y^{(p)}{\Vert }_{L^2}\to 0,\delta \to 0^{+}.$$

$$\Vert {A_n^{(p)}}^{†}{P}_n{y}^{\delta }-y^{(p)}{\Vert }_{L^2(0,2\pi )}\le C^{(p)}{n}^p\delta +({\gamma }^{(p)}+1)\Vert (I-P_n)y^{(p)}{\Vert }_{L^2},$$

$$n^{(p)}=n^{(p)}(\delta )=\kappa {\delta }^{a-\frac{1}{p}},$$

$$y^{(p)}\in H_{per}^l(0,2\pi ),l>0(\text{periodic smoothness})$$

$$\Vert {A_n^{(p)}}^{†}{P}_n{y}^{\delta }-y^{(p)}{\Vert }_{L^2(0,2\pi )}\le C^{(p)}{n}^p\delta +({\gamma }^{(p)}+1)\frac{1}{n^l}\Vert y^{(p)}{\Vert }_{H_{per}^l},$$

$$y^{(p)}\in H_{per}^{2\mu p}(0,2\pi )\mu =\frac{1}{2}\text{or}1,$$

$$y^{(p)}=a_0+\sum _{k=1}^{N_1}{b}_k\cos kx+\sum _{k=1}^{N_2}{c}_k\sin kx\in H_{per}^{\infty }(0,2\pi )$$

$$\Vert {A_n^{(p)}}^{†}{P}_{n}y{\Vert }_{L^2}⟶\infty (n\to \infty ).$$

$$\Vert {A_{n^{(p)}(\delta )}^{(p)}}^{†}{P}_{n^{(p)}(\delta )}{y}^{\delta }-y^{(p)}{\Vert }_{L^2}\to 0,\delta \to 0^{+}.$$

$$\bar{y}=y(x)-\sum _{k=0}^{p-1}\frac{y^{(k)}(0)}{k!}{x}^k\in {\mathcal{H}}_0^p(0,2\pi ),$$

$$\Vert A{x}^{†}-y\Vert =\inf \{\Vert Ax-y\Vert |x\in X\}.$$

$$A{A}^{†}y=P_{\overline{\mathcal{R}(A)}}y,\forall y\in \mathcal{R}(A)\oplus \mathcal{R}(A{)}^{\perp }.$$

$$H^p(0,2\pi ):=\{y\in L^2(0,2\pi ):D^{1}y,\cdots ,D^{p}y\in L^2(0,2\pi )\},$$

$${\int }_0^{2\pi }\zeta \phi dx=(-1{)}^k{\int }_0^{2\pi }y{\phi }^{(k)}dx,\phi \in C_0^{\infty }(0,2\pi ).$$

$$H^p(0,2\pi )={\mathcal{U}}^p[0,2\pi ]:=\{y\in C^{p-1}[0,2\pi ]:$$

$$\text{there exists}\Psi \in L^2(0,2\pi )\text{ such that}{y}^{(p-1)}(x)=\alpha +{\int }_0^x\Psi (t)dt,\alpha \in \mathrm{R}\}$$

$${\mathcal{H}}_z^p(0,2\pi ):=\{y\in H^p(0,2\pi ):y(z)=\cdots =y^{(p-1)}(z)=0\},z=0\text{or}2\pi$$

$${\dot{\mathcal{H}}}^{2p}(0,2\pi ):=\{y\in H^{2p}(0,2\pi ):$$

$$y(2\pi )=\cdots =y^{(p-1)}(2\pi )=y^{(p)}(0)=\cdots =y^{(2p-1)}(0)=0\}.$$

$$H_{per}^s(0,2\pi ):=\{\phi \in L^2(0,2\pi ):{\xi }_0^2+\sum _{k=1}^{\infty }(1+k^2{)}^s({\xi }_k^2+{\eta }_k^2)<\infty \},$$

$${\xi }_0=\frac{1}{\sqrt{2\pi }}{\int }_0^{2\pi }\phi (t)dt,{\xi }_k=\frac{1}{\sqrt{\pi }}{\int }_0^{2\pi }\phi (t)\cos ktdt,{\eta }_k=\frac{1}{\sqrt{\pi }}{\int }_0^{2\pi }\phi (t)\sin ktdt.$$

$$(\phi ,\psi {)}_{H_{per}^s}={\xi }_0{\zeta }_0+\sum _{k=1}^{\infty }(1+k^2{)}^s({\xi }_k{\zeta }_k+{\eta }_k{\lambda }_k)$$

Full text

[TABLE]

In particular, in numerical differentiations for functions above, good filtering effect (error approaches 0) is displayed with corresponding parameter choice. In addition, several numerical examples are given to show that even derivatives with discontinuity can be recovered well.

• August 2017

1 Introduction

Numerical differentiation is a classical ill-posed problem which arises in different practical fields, such as option pricing, thermodynamics and photoelectric response (See e.g. [4,11,13-16,25]). In process of numerical differentiation on a given function $y(x)$ of specific smoothness, always there would interfuse with a noise $\delta y$ in measurement or calculations. For this sake, it is routine to do numerical differentiation on the noisy function $y^{\delta}:=y+\delta y$, where the high frequency part in $\delta y$ would bring uncontrolled huge error when computing with traditional numerical algorithms. In order to overcome the difficulties of ill-posedness, several kinds of regularization method were introduced.

Tikhonov method (See [8,11,12,17,26-28]) is a classical regularization method for numerical differentiation of first order. It is generally modeled that, for $x\in[0,1]$ with a grid $\Delta=\{0=x_{0}<x_{1}<\cdots<x_{n}=1\}$ and $h=\max_{i}x_{i+1}-x_{i}$ being its mesh size, given finite noisy samples $\tilde{y}_{i}$ of $y(x_{i})$ such that $|\tilde{y}_{i}-y(x_{i})|\leq\delta$. Assume that $\tilde{y}_{0},\ \tilde{y}_{n}$ are exactly known boundary data, that is, $\tilde{y}_{0}=y(0),\ \tilde{y}_{n}=y(1)$. Then minimizing the cost functional

[TABLE]

in $\{f\in H^{2}(0,1):f(0)=y(0),f(1)=y(0)\}$ gives minimizer $f_{\alpha}$. Afterward differentiating this minimizer gives $f^{\prime}_{\alpha}$ as the regularized solution to the exact derivative $y^{\prime}$ with appropriate parameter choice $\alpha=\alpha(\delta)$. Further results illustrate that $f^{\prime}_{\alpha(\delta)}$ can converge to $y^{\prime}$ with best rate $O(\sqrt{\delta})$ with parameter choice for $\alpha=\delta^{2}$ (See [28]). However, we note that the penalty term $\|f^{\prime\prime}\|$ in cost functional basically demand that the all candidate solutions $f^{\prime}_{\alpha}$ must be at least $H^{1}$ smooth and further result in [27,28] illustrates that, for $y\in C[0,1]/H^{2}(0,1)$, under specific parameter choice $\alpha=\delta^{2}$, the upper bound for $\|f^{\prime}_{\delta^{2}}-y^{\prime}\|$ must tend to infinity as $\delta,h\to 0$. Thus this algorithm naturally deny to recover derivative with regularity less than $H^{1}$, especially discontinuous derivative.

Difference method [4,23] is another classical regularization method for numerical differentiation (including higher orders). It constructs difference quotient of $p$ order as regularized solution to exact derivatives $y^{(p)}$ with the stepsize $h$ being the regularization parameter. The convergence of this scheme is established in $L^{\infty}$ setting and will basically demand that $y^{\prime}\in C^{0,\alpha},\alpha>0$ (See [4]) which also deny to recover derivatives that are only continuous and discontinuous. Furthermore, the best convergence rate $O(\delta^{\frac{2}{3}})$ and $O(\delta^{\frac{1}{3}})$ for first and second order numerical differentiation are derived with $h=O(\delta^{\frac{1}{3}})$ respectively. But we need note the essential flaw in this algorithm that the numerical derivatives constructed by this algorithm will lose its smoothness and all be piecewise constant, whether the original function is smooth or not.

In this paper, we first formulate the $p$ order derivative $y^{(p)}$ as the unique solution of Fredholm integral equation of the first kind

[TABLE]

where $A^{(p)}:L^{2}(0,2\pi)\to L^{2}(0,2\pi)$. For the simple design of computational scheme, we apply Galerkin method with trigonometric basis to above equation to construct a regularized solutions to $y^{(p)}$ (refer to [9,19-22] for similar techniques). The basic setting can be described as below:

Assume that

[TABLE]

where $\delta$ is noise level. Given a projection sequence $\{P_{n}\}$ which project $L^{2}(0,2\pi)$ onto subspace

[TABLE]

we discrete (1.1) into a finite-rank approximation system

[TABLE]

where $A^{(p)}_{n}:=P_{n}A^{(p)}P_{n}:X_{n}\longrightarrow X_{n}$.

Now solving (1.2) in sense of Moore-Penrose inverse gives ${A_{n}^{(p)}}^{\dagger}P_{n}y$. This is a natural approximation scheme for $y^{(p)}$, where $\dagger$ denotes the Moore-Penrose inverse of linear operator. Considering the noise $\delta y$, above scheme should be adjusted to a regularized scheme

[TABLE]

where $n^{(p)}:=n^{(p)}(\delta)$ is the regularization parameter choice such that

[TABLE]

and

[TABLE]

Here notice that $n^{(p)}:=n^{(p)}(\delta)\ (p=1,2,3)$ stands for parameter choice strategy of the first three order numerical differentiation respectively. Throughout this paper, without special indication, we follow this notation and $p\in\overline{1,2,3}$.

Main results of this paper and corresponding remarks are listed as follows.

• •

For $y\in{\mathcal{H}}^{p}_{0}(0,2\pi)$ (this restriction on initial value data is removable), where ${\mathcal{H}}^{p}_{0}(0,2\pi):=\{y\in H^{p}(0,2\pi):y(0)=\cdots=y^{(p-1)}(0)=0\}.$ a priori error estimate is obtained uniformly for first three order numerical differentiation as

[TABLE]

This determines the parameter choice strategy:

[TABLE]

where $a\in(0,\frac{1}{p})$ is optional and $\kappa$ is a constant which depends on the concrete form of $\|(I-P_{n})y^{(p)}\|_{L^{2}}$. This establish a convergence result for numerical differentiation of first three order when $y^{(p)}\in L^{2}$, especially for derivative with discontinuities. However, we need specify that, when recovering $y^{(p)}\in L^{2}$ are only continuous and discontinuous with no periodic smoothness, the constant $\kappa$ is unknown and need to test in experiments (See section 8.3). In addition we give a notice that, whether the derivative is smooth or not, its approximation by above algorithm will be real analytic since it is a trigonometric polynomial.

• •

Supplemented with a priori information

[TABLE]

above error estimate is strengthened into a more explicit form as

[TABLE]

where $C^{(p)}$, $\gamma^{(p)}$ are all independent constants given in proceeding sections. And optimal convergence rate $O(\delta^{\frac{2\mu}{2\mu+1}})$ can be derived under periodic Sobolev source condition

[TABLE]

with parameter choice $n^{(p)}=\lambda^{(p)}\delta^{-\frac{1}{(2\mu+1)p}}$, where $\lambda^{(p)}$ is a constant only depends on exact derivative $y^{(p)}$ and can be given explicitly in preceding sections 6,7. In particular, when

[TABLE]

the optimal parameter choice will degenerate to a constant $n=\max(N_{1},N_{2})$ which does not depend on noise. Furthermore, the numerical study in section 8.1 demonstrates the good filtering effect (error approaches 0) occurs in this specific case.

• •

In a more general setting for $p$ order numerical differentiation when $y\in H^{p}(0,2\pi)$, it is indicated in Corollary 6.2 that, when $y\in H^{p}(0,2\pi)\setminus{\mathcal{H}}^{p}_{0}(0,2\pi)$,

[TABLE]

Now any parameter choice $n^{(p)}(\delta)$ such that $n^{(p)}=n^{(p)}(\delta)\to\infty\ (\delta\to 0^{+})$ may not be a proper regularization parameter since we can not determine

[TABLE]

through traditional estimate any more (See second point behind Corollary 5.2). In order to recover the regularization effect of algorithm, we introduce Taylor polynomial truncation of $p-1$ order to reform the regularized scheme, that is, using

[TABLE]

to replace $y\in H^{p}(0,2\pi)$. In this way, the regularization effect can be well recovered (See section 6,7)with exact measurements on initial value data. Furthermore, we take possible noise in measurements in initial value data into consideration, and this effectively relax the requirement on precision of initial value data.

Outline of Paper: In section 2, we introduce some tools and basic lemmas. In section 3, we illustrate general framework, give the main idea on how to utilize the noisy data $y^{\delta}$ to recover the $p$ order derivatives $y^{(p)}$. In section 4, we give corresponding analytic solution formula to Galerkin approximation system which determines the well-posedness result and upper bound for noise error. In section 5, we propose an estimate on approximation error when RHS $y$ belongs to $\mathcal{H}^{p}_{0}(0,2\pi)$, and give the convergence and divergence results with respect to $y\in\mathcal{H}^{p}_{0}(0,2\pi)$ and $y\in L^{2}(0,2\pi)\setminus\mathcal{H}^{p}_{0}(0,2\pi)$ respectively. In sections 6 and 7, with periodic Sobolev source condition of order $2\mu p$, we construct a priori error estimate and indicate the parameter choice strategy for optimal convergence rate $O(\delta^{\frac{2\mu}{2\mu+1}})$ when $y\in\mathcal{H}^{p}_{0}(0,2\pi)$ and $y\in H^{p}(0,2\pi)\setminus\mathcal{H}^{p}_{0}(0,2\pi)$ respectively. In section 8, we test some numerical examples to show the characteristics and effects of algorithm when derivatives are smooth and discontinuous respectively. In section 9, we conclude the main work of this paper.

2 Preliminary and Basic Lemmas

2.1 Moore-Penrose inverse

Let $X,Y$ be Hilbert space, and $A$ be bounded linear operator mapping from $X$ to $Y$. $\mathcal{D}(A)$, $\mathcal{N}(A)$ and $\mathcal{R}(A)$ denote its domain, null space and range, respectively.

For $A:X\to Y$ and $y\in\mathcal{R}(A)\oplus\mathcal{R}(A)^{\perp}$, the Moore-Penrose inverse $x^{\dagger}:=A^{\dagger}y$ is defined as the element of smallest norm satisfying

[TABLE]

Thus $A^{\dagger}:\mathcal{D}(A^{\dagger}):=\mathcal{R}(A)\oplus\mathcal{R}(A)^{\perp}\subseteq Y\longrightarrow X$ defines a closed linear operator from $Y$ to $X$.

In the following, we indicate some useful properties of Moore-Penrose inverse $A^{\dagger}$:

• •

If $A:X\to Y$ is one-to-one, then, for $y\in\mathcal{R}(A)$, $A^{\dagger}y$ naturally degenerates into $A^{-1}y$.

• •

If $\mathcal{R}(A)$ is closed, then $\mathcal{D}(A^{\dagger})=\mathcal{R}(A)\oplus\mathcal{R}(A)^{\perp}=Y$ and by closed graph theorem, $A^{\dagger}:Y\to X$ is bounded.

• •

If $\mathcal{R}(A)$ is closed, then $AA^{\dagger}=P_{\mathcal{R}(A)},\ A^{\dagger}A=P_{\mathcal{N}(A)^{\perp}}.$ If $\mathcal{R}(A)$ is not necessarily closed, then the former identity need be adjusted into

[TABLE]

For more comprehensive information on Moore-Penrose inverses, see [2, Chapter 9] or [6,7].

2.2 Sobolev spaces

Throughout this paper, we only discuss on Sobolev space over $\mathrm{R}$. Without specification, we denote $H^{p}(0,2\pi):=H_{\mathrm{R}}^{p}(0,2\pi)$. Here we introduce all kinds of notations of Sobolev spaces which will be used in the context. For more information, see [1,5] and [3, Appendix 4].

2.2.1 Sobolev spaces of integer order

For some positive integer $p$, the Sobolev space $H^{p}(0,2\pi)$ is defined as

[TABLE]

where $D^{k}y$ means weak derivative, defined as $\zeta\in L^{2}(0,2\pi)$ which satisfies

[TABLE]

Equivalently, it can be characterized in absolute continuous form (refer to [3, Page 14]) as

[TABLE]

[TABLE]

Here notice that above $"="$ admits a possible change in a set of measure zero. In this paper, when it concerns Sobolev functions of one variable $y\in H^{p}(0,2\pi)$, we, by default, modify $y\in H^{p}(0,2\pi)\ (p\in\mathrm{N})$ in a set of measure zero such that it belongs to the latter fine function space $\mathcal{U}^{p}[0,2\pi]$.

Besides, for $p\in\mathrm{N}$, we define

[TABLE]

and

[TABLE]

[TABLE]

2.2.2 Fractional periodic Sobolev spaces

For real number $s>0$, periodic Sobolev spaces of fractional order $H_{per}^{s}(0,2\pi)$ is defined in trigonometric form as

[TABLE]

where

[TABLE]

Supplementing another element $\psi\in H_{per}^{s}(0,2\pi)$, its inner product is rephrased as

[TABLE]

with

[TABLE]

In addition, we define

[TABLE]

2.3 Integro-differential operator of $p$ order

Define integro-differential operator of integer order $p$ as:

[TABLE]

[TABLE]

This is a compact linear operator with infinite-dimensional range, which satisfies

Lemma 2.1

[TABLE]

Proof 1

”$\subseteq$”: Assume that $y\in{\mathcal{H}}^{p}_{0}(0,2\pi)\ (p=1,2,3)$. There exists a $y^{\star}\in\mathcal{U}^{p}[0,2\pi]$ as a modification of $y$ in a possible [math] measure set with [math] initial value data, that is, $y^{\star}(0)=\cdots={y^{\star}}^{(p-1)}(0)=0$. With integration formula by parts, it is not difficult to verify that,

[TABLE]

Thus,

[TABLE]

”$\supseteq$”: For simplicity, we only provide proof of case $p=2$. Assuming $y\in\mathcal{R}(A^{(2)})$, then there exists a $\varphi\in L^{2}(0,2\pi)$ such that

[TABLE]

It is not difficult to verify that

[TABLE]

With definition of (2.2), it yields that $y\in H^{2}(0,2\pi)$. Then by absolute continuous characterization (2.3) of Sobolev function of one variable, there exist a $y^{\star}\in\mathcal{U}^{2}[0,2\pi]$ as modification of $y$ in a [math] measure set. Thus we have

[TABLE]

Notice that

[TABLE]

are all continuous functions, thus

[TABLE]

”$\supseteq$” holds for case $p=2$.

With above equality, we describe the density of range in $L^{2}(0,2\pi)$.

Lemma 2.2

[TABLE]

where $I$ is the identity operator on $L^{2}(0,2\pi)$.

Proof 2

With Lemma 2.1, ${\mathcal{H}}^{p}_{0}(0,2\pi)=\mathcal{R}(A^{(p)}).$ Recall the fact that $C^{\infty}_{0}(0,2\pi)$ is dense in $L^{2}(0,2\pi)$ and notice that $C^{\infty}_{0}(0,2\pi)\subseteq{\mathcal{H}}^{p}_{0}(0,2\pi)$, then $L^{2}(0,2\pi)=\overline{C^{\infty}_{0}(0,2\pi)}\subseteq\overline{{\mathcal{H}}^{p}_{0}(0,2\pi)}\subseteq L^{2}(0,2\pi).$ The result follows.

This implies that $\mathcal{R}(A^{(p)})^{\perp}={0}$, and

[TABLE]

Now differentiating the both sides of equation (1.1) for $y\in\mathcal{H}^{p}_{0}(0,2\pi)$ in $p$ order yields that ${A^{(p)}}^{\dagger}y={A^{(p)}}^{-1}y=y^{(p)}$. This gives

Lemma 2.3

${A^{(p)}}^{\dagger}y=y^{(p)},\ \forall y\in{\mathcal{H}}^{p}_{0}(0,2\pi).$**

2.4 Galerkin Projection scheme with Moore-Penrose inverses

Let $X$ be Hilbert space. For the linear operator equation

[TABLE]

where $A:X\longrightarrow X$ is bounded and linear. To approximate

[TABLE]

We introduce a sequence of finite-dimensional subspaces $\{X_{n}\}$, which satisfies

[TABLE]

Then construct a sequence of orthogonal projections $\{P_{n}\}$, where $P_{n}$ projects $X$ onto $X_{n}$, and gives Galerkin approximation setting

[TABLE]

where $A_{n}:=P_{n}AP_{n}:X_{n}\longrightarrow X_{n}$. Hence solving (2.5) in sense of Moore-Penrose inverse gives Galerkin projection scheme

[TABLE]

where $A^{\dagger}_{n}:\mathcal{R}(A_{n})+{\mathcal{R}(A_{n})}^{\perp_{n}}=X_{n}\longrightarrow X_{n}$. Notice that $\perp_{n}$ means orthogonal complement in finite dimensional Hilbert space $X_{n}$.

Now $\{\varphi^{\dagger}_{n}\}$ is a natural approximate scheme for $\varphi^{\dagger}$. To study its convergence property, we introduce the Groetsch regularizer for setting (2.5) as

[TABLE]

define the Groetsch regularity as $\ \sup\limits_{n}\|R_{n}\|<+\infty$, and introduce the following result:

Lemma 2.4

*For above Galerkin approximate setting (2.5), if Groetsch regularity holds, then

(a) For $y\in\mathcal{D}(A^{\dagger})=\mathcal{R}(A)+\mathcal{R}(A)^{\perp}$.*

[TABLE]

(b) For $y\notin\mathcal{D}(A^{\dagger})$,

[TABLE]

Proof 3

see[18, theorem 2.2]

2.5 Higher order estimate under trigonometric basis

To further estimate the right term $\|(I-P_{n})A^{\dagger}y\|$ under trigonometric basis in $L^{2}$, similar to the result [3, Lemma A.43], we introduce another error estimate:

Lemma 2.5

Let $P_{n}:L^{2}(0,2\pi)\longrightarrow X_{n}\subset L^{2}(0,2\pi)$ be an orthogonal projection operator, where

[TABLE]

Then $P_{n}$ is given as follows

[TABLE]

where

[TABLE]

[TABLE]

are the Fourier coefficients of $x$. Furthermore, the following estimate holds:

[TABLE]

where $r\geq 0$.

3 General Framework

We start from

Problem 3.1

Assume that we have $y\in{\mathcal{H}}^{p}_{0}(0,2\pi)$ and $y^{\delta}$ measured on $(0,2\pi)$, belonging to $L^{2}(0,2\pi)$ such that $\|y^{\delta}-y\|_{L^{2}}\leq\delta$. How to get a stable approximation to $y^{(p)}$ ?

In Lemma 2.3, we have known that $y^{(p)}$ is the solution of linear operator equation

[TABLE]

when

[TABLE]

3.1 Formulation of finite-dimensional approximation system

In the following, we consider to approximate $y^{(p)}$ by the Galerkin method. Set

[TABLE]

Choose a sequence of orthogonal projection operators $\{P_{n}\}$, where $P_{n}$ projects $L^{2}(0,2\pi)$ onto

[TABLE]

Then degenerate the original operator equation with noisy data

[TABLE]

into a finite-rank system

[TABLE]

where

[TABLE]

Span $A_{n}^{(p)}$ under above basis, then the finite-rank system (3.2) is transformed into the linear system as

[TABLE]

Notice that $M_{n}^{(p)}$ and $b^{\delta}_{n}$ are defined as follows:

[TABLE]

where

[TABLE]

[TABLE]

Indeed,

[TABLE]

[TABLE]

And $b^{\delta}_{n}:=(f_{0},f_{1},g_{1},\cdots,f_{n},g_{n})^{T}$ is defined as

[TABLE]

[TABLE]

Indeed,

[TABLE]

Once we figure out ${u_{n}^{p,\delta}}={M_{n}^{(p)}}^{\dagger}b^{\delta}_{n}$, then we obtain solution for (3.2),

[TABLE]

in sense of Moore-Penrose inverse, ${\varphi_{n}^{p,\delta}}={A_{n}^{(p)}}^{\dagger}y^{\delta}_{n}$. This is the regularized scheme. In the following, we need to determine a regularization parameter $n^{(p)}=n^{(p)}(\delta)$ such that

[TABLE]

3.2 Total error estimate and parameter choice for regularization

Now, in order to control the accuracy of computation, we adjust parameter choice strategy $n^{(p)}=n^{(p)}(\delta)$ according to following total error estimate

[TABLE]

Since Lemma 2.3 illustrates that

[TABLE]

inserting (3.7) into (3.6), the formula (3.6) becomes

[TABLE]

Throughout this paper we use the following definitions

• •

Total error

[TABLE]

which is broken into two parts (c.f.[10, Chapter 1.1]):

• •

Noise error:

[TABLE]

• •

Approximation error:

[TABLE]

It is an easy observation that $e^{(p)}_{T}\leq e^{(p)}_{N}+e^{(p)}_{A}$. Upon this fact, we figure out the total error estimate by estimating $e^{(p)}_{N}$ and $e^{(p)}_{A}$ respectively.

4 Well-posedness and numerical scheme of Galerkin System

With concrete expressions of $M^{(p)}_{n}$ in Appendix A, it is not difficult to obtain:

Theorem 4.1

Finite dimensional system (3.4) is well-posed, that is, there exists a unique solution to (3.4), denoted as

[TABLE]

where

[TABLE]

Moreover, analytic formulas for the solution of Galerkin approximation system (3.2) are determined as follows:

[TABLE]

*Corresponding three cases are listed as follows.

Case $p=1$:*

[TABLE]

[TABLE]

[TABLE]

Case $p=2$:

[TABLE]

[TABLE]

[TABLE]

where

[TABLE]

Case $p=3$:

[TABLE]

[TABLE]

[TABLE]

where

[TABLE]

[TABLE]

Remark 4.1

After we solve the Galerkin approximation system, we know that $A^{(p)}_{n}:X_{n}\to X_{n}$ is one-to-one and surjective. For the usage of the proceeding section, we claim that $\mathcal{N}(A^{(p)}_{n})=0,\ \mathcal{R}(A^{(p)})=X_{n}.$

Remark 4.2

With above analytic formulas, it is not difficult to figure out that

[TABLE]

where $C^{(1)}=\sqrt{3},\ C^{(2)}\approx 11.8040,\ C^{(3)}\approx 345.0754$. Here we specify that $\|\cdot\|_{X_{n}}$ is induced by $\|\cdot\|_{L^{2}}$, that is, $\|x_{n}\|_{X_{n}}:=\|x_{n}\|_{L^{2}}$, $\forall x_{n}\in X_{n}$. This give bound to the estimate of noise error.

5 Estimate on Approximation Error and Instability result

We use Lemma 2.4 to analyse the convergence and divergence of Galerkin method. The key point is the estimate of

[TABLE]

To gain an uniform upper bound for above formula, we first prepare two decay estimate of

[TABLE]

with respect to integer variable $j$:

Lemma 5.1

For operators $A^{(p)},A^{(p)}_{n}$ defined in (2.4),(3.3) respectively, set

[TABLE]

When $j\geq n+1$,

[TABLE]

[TABLE]

[TABLE]

and

[TABLE]

where

[TABLE]

When $p=3$, we need an extra condition $n\geq 5$ to maintain above estimate.

Proof 4

Case $p=1$: When $j\geq n+1$, substituting (B.1) into (4.1),(4.2) and (4.3), it follows that

[TABLE]

*This gives lemma for case $p=1$.

** Case $p=2$**: Inserting (B.2) into (4.4), it follows that*

[TABLE]

Hence

[TABLE]

Besides, inserting (B.2) into (4.5),(4.6) respectively, it yields that

[TABLE]

Then, by (5.3),

[TABLE]

Case $p=3$: Inserting (B.3) into (4.7), it follows that

[TABLE]

Notice Proposition C.1 (C.4),

[TABLE]

Hence,

[TABLE]

Besides, insert (B.3) into (4.8), then it follows that

[TABLE]

By Proposition C.1(C.2), it is routine to obtain that

[TABLE]

Hence, with (5.4),

[TABLE]

Further, insert (B.3) into (4.9), and we have

[TABLE]

Hence

[TABLE]

Lemma 5.2

For operators $A^{(p)},A^{(p)}_{n}$ defined in (2.4),(3.3) respectively, set

[TABLE]

When $j\geq n+1$,

[TABLE]

[TABLE]

[TABLE]

and

[TABLE]

where

[TABLE]

Notice that when $p=3$, we need the extra condition $n\geq 5$ to maintain above estimate.

Proof 5

Case $p=1$: When $j\geq n+1$, insert (B.4) into (4.1),(4.2),(4.3), then

[TABLE]

*This gives lemma for case $p=1$. *

Case $p=2$: Insert (B.5) into (4.4), and it follows that

[TABLE]

With Proposition C.1 (C.1), it follows that

[TABLE]

Besides, insert (B.5) into (4.5),(4.6), then

[TABLE]

Then by (5.7) we have

[TABLE]

Case $p=3$: Insert (B.6) into (4.7), and it follows that

[TABLE]

Notice that it is easy to obtain that

[TABLE]

from Proposition C.1 (C.3). In this way, with Proposition C.1 (C.4), when $n\geq 5$,

[TABLE]

[TABLE]

Besides, insert (B.6) into (4.8), and we have

[TABLE]

Hence, by (5.8)

[TABLE]

Further, insert (B.6) into (4.9), and it follows that

[TABLE]

Hence, by (5.8)

[TABLE]

Lemma 5.3

Set $A^{(p)},A_{n}^{(p)}$ defined in (2.4), (3.3) respectively. Then

[TABLE]

where

[TABLE]

Remark 5.1

With direct computations, we can obtain that

[TABLE]

Proof 6

*Set $v=a_{0}\frac{1}{\sqrt{2\pi}}+\sum^{\infty}_{k=1}a_{k}\frac{\cos kt}{\sqrt{\pi}}+\sum^{\infty}_{k=1}b_{k}\frac{\sin kt}{\sqrt{\pi}}$ such that $\|v\|_{L^{2}}=1$, that is, $a^{2}_{0}+\sum^{\infty}_{k=1}a^{2}_{k}+\sum^{\infty}_{k=1}b^{2}_{k}=1.$ We consider the estimate on $\|K^{(p)}_{n}v\|_{L^{2}},\ n\in\mathrm{N}$.

Since ${A^{(p)}_{n}}^{\dagger}P_{n}A^{(p)}$ $(n\in\mathrm{N})$ is continuous with Remark 4.2,*

[TABLE]

[TABLE]

[TABLE]

Recall Lemma 5.1 and Lemma 5.2. It follows that

[TABLE]

where

[TABLE]

[TABLE]

[TABLE]

By (5.2), (5.6) and the Cauchy inequality, we have

[TABLE]

[TABLE]

[TABLE]

(5.10),(5.11),(5.12) together with (5.9) give that, for all $v$ such that $\|v\|_{L^{2}}=1$,

[TABLE]

[TABLE]

where $C^{(p)}_{i}\ (p=1,2,3;\ i=1,2,\cdots,6)$ are all constants defined in Lemma 5.1 and Lemma 5.2.

Theorem 5.1

The Groetsch regularity holds for Galerkin setting (3.2); that is,

[TABLE]

where

[TABLE]

and $A^{(p)},A_{n}^{(p)}$ are defined in (2.4), (3.3) respectively.

Proof 7

Since

[TABLE]

[TABLE]

by Lemma 5.3 we have

[TABLE]

After the examination of (5.1), we have an estimate on the approximation error.

Corollary 5.1

For ${A^{(p)}_{n}}$ defined as (3.3), we have

[TABLE]

for every $y\in{\mathcal{H}}^{p}_{0}(0,2\pi)$. Furthermore, with a priori information $y^{(p)}\in H_{per}^{l}(0,2\pi)$, it yields that

[TABLE]

where $\gamma^{(p)}$ is constant given in Theorem 5.1.

Proof 8

By Lemma 2.2, 2,3, for $y\in{\mathcal{H}}^{p}_{0}(0,2\pi)$,

[TABLE]

Using Lemma 2.4, Remark 4.1, Theorem 5.1, it yields that

[TABLE]

Now provided with a priori information $y^{(p)}\in H_{per}^{l}(0,2\pi),\ l>0.$ It yields that

[TABLE]

from Lemma 2.5. This improves (5.13) to the latter result needed.

Corollary 5.2

For ${A^{(p)}_{n}}$ defined as (3.3), we have

[TABLE]

for every $y\in L^{2}(0,2\pi)\setminus\mathcal{R}(A^{(p)})$, where $\mathcal{R}(A^{(p)})={\mathcal{H}}^{p}_{0}(0,2\pi)$.

Proof 9

Lemma 2.4 tells that $\mathcal{D}({A^{(p)}}^{\dagger})=\mathcal{R}(A^{(p)})$. Since the estimate of (5.1) holds, with Lemma 2.4 (b), the result surely holds.

Here Corollary 5.2 tells us two questions:

• •

the first question is, for $p$ order numerical differentiation, when $y\in{\mathcal{H}}^{p}_{0}(0,2\pi)$, with interfuse of noise $\delta y$, $y^{\delta}$ would generally locate in $L^{2}(0,2\pi)\setminus\mathcal{H}^{p}_{0}(0,2\pi)$. Then with increasing choice of index $n$ independent of noise level $\delta$,

[TABLE]

This fact shows that, without proper parameter choice strategy for $n^{(p)}:=n^{(p)}(\delta)$, numerical scheme constructed as ${A^{(p)}_{n}}^{\dagger}P_{n}y^{\delta}$ is natively instable.

• •

the second question is a worse one. With the more general setting $y\in H^{p}(0,2\pi)$ for $p$ order Numerical differentiation, if $y\in H^{p}(0,2\pi)\setminus{\mathcal{H}}^{p}_{0}(0,2\pi)$, then, with any parameter choice strategy $n^{(p)}:=n^{(p)}(\delta)$ such that $n^{(p)}(\delta)\to+\infty\ (\delta\to 0^{+})$, approximation error

[TABLE]

In addition with estimate on noise error $e^{(p)}_{N}\leq C^{(p)}n^{(p)}(\delta)\to\infty\ (\delta\to 0^{+})$ (by (4.10)), one could see the invalidness of single regularization parameter choice since only adjusting parameter choice $n=n^{(p)}(\delta)$ can not gives $e^{(p)}_{T}\to 0\ (\delta\to 0^{+})$.

The following two sections will answer above two questions respectively.

6 Total Error Estimate for $y\in{\mathcal{H}}^{p}_{0}(0,2\pi)$ and Parameter Choice for Regularization

To solve the first question, we introduce regularization in the following procedure:

Combining Corollary 5.1 with Remark 4.2, it gives

Theorem 6.1

Set ${A^{(p)}_{n}}$ as (3.3), $y\in{\mathcal{H}}^{p}_{0}(0,2\pi)$. Then

[TABLE]

Furthermore, with a priori information $y^{(p)}\in H_{per}^{l}(0,2\pi)$,

[TABLE]

Remark 6.1

In the case that $y\in{\mathcal{H}}^{p}_{0}(0,2\pi)$ is provided but no a priori information on exact solution $y^{(p)}$, we determine parameter choice strategy from (6.1) as

[TABLE]

*where $a\in(0,\frac{1}{p})$ is optional. However, we specify that, in this case, the convergence rate can not be obtained higher than $O(1)$.

In the case that $y\in{\mathcal{H}}^{p}_{0}(0,2\pi)$ and $y^{(p)}\in H^{l}_{per}(0,2\pi)$, we could determine parameter choice strategy from (6.2) as*

[TABLE]

Hence it follows that

[TABLE]

where

[TABLE]

Remark 6.2

Assume that

[TABLE]

Choosing $n_{3}^{(p)}=n_{3}^{(p)}(\delta)=\delta^{-\frac{1}{(2\mu+1)p}}$, we gain the optimal convergence rate from (6.5), that is,

[TABLE]

Here we specify that (6.6) is a slightly variant version of the standard source condition stated in [10], that is, $y^{(p)}\in\mathcal{R}({A^{(p)}}^{*}A^{(p)})^{\mu}\ (\mu=\frac{1}{2}\ \textrm{or}\ 1)$, where

[TABLE]

Notice that

[TABLE]

and

[TABLE]

Remark 6.3

Assume that noise level $\delta$ range in any closed interval $[\delta_{0},\delta_{1}]\subseteq(0,+\infty)$. If $y\in{\mathcal{H}}^{p}_{0}(0,2\pi)$ and $y^{(p)}\in H_{per}^{\infty}(0,2\pi)$ such that $\lim_{l\to\infty}\|y^{(p)}\|^{\frac{1}{l+p}}_{H_{per}^{l}}$ exists, then the regularization parameter is determined by (6.4) as

[TABLE]

which only depends on exact derivative $y^{(p)}$, not concerned with noise level $\delta$. Furthermore, in the case that $y\in\mathcal{H}^{p}_{0}(0,2\pi)$ and

[TABLE]

where $b_{N_{1}},\ c_{N_{2}}\neq 0$. The regularization parameter is determined by (6.7) as

[TABLE]

which only depends on the highest frequency of trigonometric polynomial in (6.8), not concerned with noise level $\delta$.

7 Extended Numerical Differentiation on $H^{p}(0,2\pi)$

7.1 Extended result with exact measurements at endpoint $x=0$

Theorem 6.1 provides a result of stable numerical differentiation on $y\in{\mathcal{H}}^{p}_{0}(0,2\pi)$, where

[TABLE]

We consider to remove the restriction on initial value data, and extend the result into case $y\in H^{p}(0,2\pi)$.

Observing that, for $y\in H^{p}(0,2\pi)$,

[TABLE]

we naturally adjust regularized scheme (1.3) into

[TABLE]

Now, given exact measurements on the initial value data,

[TABLE]

we can adjust Theorem 6.1 into the following version.

Theorem 7.1

Set ${A^{(p)}_{n}}$ as (3.3), $y\in H^{p}(0,2\pi)$. Then

[TABLE]

Furthermore, with a priori information $y^{(p)}\in H_{per}^{l}(0,2\pi)$,

[TABLE]

7.2 Extended result with noisy measurements at endpoint $x=0$

However, in practical cases, one can not obtain initial value data $y(0),y^{\prime}(0),y^{\prime\prime}(0)$ exactly. Instead, one could only obtain a cluster of noisy data, denoted as $\Lambda_{0}(0),\Lambda_{1}(0),\Lambda_{2}(0)$ respectively. Now provided with above endpoint measurement, we reformulate the problem of $p$ order numerical differentiation as:

Problem 7.1

Assume that we have

• •

$y\in H^{p}(0,2\pi)$* and $y^{\delta}$ measured on $(0,2\pi)$, which belongs to $L^{2}(0,2\pi)$ such that $\|y^{\delta}-y\|_{L^{2}}\leq\delta$,*

• •

Noisy initial value data $\Lambda_{0}(0),\ \cdots,\Lambda_{p-1}(0)$ for $y(0),\cdots,y^{(p-1)}(0)$ respectively, which satisfies that

[TABLE]

How to gain stable approximation to $y^{(p)}$?

An estimate similar to (6.2) is constructed to answer this question:

Theorem 7.2

Set ${A^{(p)}_{n}}$ as (3.3), $y\in H^{p}(0,2\pi)$ and $y^{(p)}\in H_{per}^{l}(0,2\pi)$. Then

[TABLE]

[TABLE]

For convenience of notations, we set $\delta_{i}=\delta$, then it follows that

[TABLE]

[TABLE]

where $C_{\Delta}^{(p)}:=(\Delta_{p}+1)C^{(p)}$ and $\Delta_{p}:=\sum\limits^{p-1}_{k=0}\frac{\|x^{k}\|_{L^{2}}}{k!}$.

Remark 7.1

In this case, it is necessary to specify that the parameter choice strategy should be adjusted from (6.4) to the following,

[TABLE]

Also, assume that noise level $\delta$ range in the interval $(\delta_{0},\delta_{1})$, where $0<\delta_{0}<<1$ and $\delta_{1}<+\infty$. When $y\in H^{p}(0,2\pi)$ and $y^{(p)}\in H_{per}^{\infty}$ such that $\lim_{l\to\infty}\|y^{(p)}\|^{\frac{1}{l+p}}_{H_{per}^{l}}$ exists, the regularization parameter is determined by (7.1) as

[TABLE]

which remains not concerned with noise level $\delta$, also not concerned with the additional noise level $\delta_{i}$ of initial value data. Besides, in the case that

[TABLE]

where $b_{N_{1}},\ c_{N_{2}}\neq 0$. The optimal parameter choice is determined by (7.2) as

[TABLE]

which is just the same as (6.9), still not concerned with noise level $\delta$ and additional noise in initial value data.

Remark 7.2

The optimal convergence rate $O(\delta^{\frac{2\mu}{2\mu+1}})$ can be achieved in the same way as Remark 6.2.

Proof 10

For $y\in H^{p}(0,2\pi)$ and $y^{\delta}\in L^{2}[0,2\pi]$ with $\|y^{\delta}-y\|\leq\delta$,

[TABLE]

where

[TABLE]

[TABLE]

Apply Theorem 7.1, and it follows that

[TABLE]

Besides,

[TABLE]

where $\Delta_{p}:=\sum\limits^{p-1}_{k=0}\frac{\|x^{k}\|_{L^{2}}}{k!}$. Then we have

[TABLE]

[TABLE]

8 Numerical Experiments

All experiments are performed in Intel(R) Core(TM) i7-7500U CPU @2.70GHZ 2.90 GHZ Matlab R 2017a. For all experiments, the regularized solution is given by

[TABLE]

with regularization parameter choice $n=n^{(p)}=n^{(p)}(\delta,\delta_{i})(p=1,2,3)$, where

[TABLE]

and

[TABLE]

All experiments are divided into two cases:

• •

Case I: $\delta\neq 0,\delta_{i}=0$, that is, high frequency noise $\delta y$ and exact initial value data.

• •

Case II: $\delta=\delta_{i}\neq 0$, that is, high frequency noise $\delta y$ and noisy initial value data.

The following index is introduced to measure the computational accuracy in tests:

• •

Relative error

[TABLE]

8.1 On smooth functions

Example 8.1

Set

[TABLE]

[TABLE]

We use the $y^{\delta}_{i}$ as test function for $i$ order numerical differentiation. Notice that

[TABLE]

[TABLE]

[TABLE]

[TABLE]

8.1.1 Unified observation on cases with smooth derivative

We first investigate into the case $(\delta,\delta_{i})=(0.01,0)$. All data in this case can be divided into three phases $\mathcal{V}_{1}=\{2,4\}$, $\mathcal{V}_{2}=\{6,8\}$, $\mathcal{V}_{3}=\{12\}$.

We can compare above three phases and quickly find that when $n\in\mathcal{V}_{2}$, the relative error $r$ is the least and almost approaches [math]. This displays the good filtering effect of algorithm on specific class of functions (the sum of trigonometric polynomial and polynomial of order less than order $p$, where $p$ denotes the order of numerical differentiation), and also correspond to the fact indicated by (6.9): the best regularization parameter $n=6$.

Now we explain the source of the good filtering effect in case $p=2$, the other cases are similar.

• •

The exact measurement of initial value data help give a precise Taylor polynomial truncation to eliminate the polynomial term $q_{2}(x)$ in $y^{\delta}_{2}$, thus the computational scheme is transformed into ${A^{(2)}_{n}}^{\dagger}P_{n}(\bar{p}(x)+\delta y)$, where $\bar{p}(x)=p(x)-\sum^{6}_{k=1}\frac{1}{k}x$.

• •

The parameter choice $n=6$ appropriately eliminate the noise component $\delta y$, now ${A^{(2)}_{6}}^{\dagger}P_{6}(\bar{p}(x)+\delta y)={A^{(2)}_{6}}^{\dagger}P_{6}\bar{p}(x)$. Notice that $\bar{p}(x)\in\mathcal{R}(A^{(2)})$, ${A^{(2)}}^{\dagger}\bar{p}(x)=y^{\prime\prime}_{2}$,

[TABLE]

[TABLE]

Then the uniqueness of Galerkin system in Theorem 4.1 gives that ${A^{(2)}_{6}}^{\dagger}P_{6}\bar{p}(x)=y^{\prime\prime}_{2}$, that is, approximate solution strictly equals to the exact solution in this case, the good filtering effect appears.

The deviation of the accuracy in $\mathcal{V}_{1}$ and $\mathcal{V}_{3}$ can also be explained in the same way. The former case of choice $n=2$ with lower accuracy is due to that $P_{n}y^{\delta}_{2},\ n=2$ does not cover the major part of $p(x)$, but with the increase of $n\in\mathcal{V}_{1}$, the coverage increases and hence the accuracy improves. As to the latter case $\mathcal{V}_{3}$, now the $\delta y=\frac{\sin 12x}{\sqrt{\pi}}$ come into the computation, thus the good filtering effect disappears.

For case II with noise pair $(\delta,\delta_{i})=(0.01,0.01)$, it can be seen in Table 1 that, when regularization parameter $n\in\mathcal{V}_{2}$, the corresponding relative error $r$ all approach or exceed 0.01 uniformly. Now, compared to the case I with the same choice for regularization parameter, the good filtering effect of case I are strongly weakened. This is because the noise $\delta_{i}$ in initial value data bring a not complete Taylor polynomial truncation and hence parameter $n=6$ can not eliminate the lower-frequency noise components in

[TABLE]

8.2 On periodic weakly differentiable derivatives

In this subsection, we mainly investigate the effectiveness of parameter choice (6.4) and (7.1) for general case (compared with the case in subsection 8.1).

8.2.1 First order

Example 8.2

Set

[TABLE]

[TABLE]

[TABLE]

[TABLE]

8.2.2 Second order

Example 8.3

Set

[TABLE]

[TABLE]

[TABLE]

[TABLE]

8.2.3 Third order

Example 8.4

Set

[TABLE]

[TABLE]

[TABLE]

[TABLE]

8.2.4 Unified observation on cases with periodic weakly differentiable derivative

In this subsection, we utilize strategies (6.4) and (7.1) to determine regularization parameter for two cases with different noise level respectively. The figures 1-3 show the good effectiveness of strategy proposed in this paper in a general aspect.

8.3 On discontinuous derivatives

In the following numerical examples with non-periodic discontinuous derivatives, we choose to adjust parameter

[TABLE]

with experiments, not by (6.3) for the uncertainties to determine $\kappa$. Figures corresponding to least-error $r$ in numerical differentiation of each order are attached.

8.3.1 First order

Example 8.5

Set

[TABLE]

[TABLE]

[TABLE]

[TABLE]

8.3.2 Second order

Example 8.6

Set

[TABLE]

[TABLE]

[TABLE]

[TABLE]

8.3.3 Third order

Example 8.7

Set

[TABLE]

[TABLE]

and

[TABLE]

[TABLE]

8.3.4 Unified observation on cases with discontinuous derivative

It can be concluded from figure 4,5,6 that, in both cases, when regularization parameters are chosen appropriately, the computational error can be well controlled. However, for sake of the intersection of frequency band of $y$ and $\delta y$ (this does not happen in the example we list in the subsection 8.1), the good filtering effect disappears in case with discontinuous derivative. Besides, we note that when the choice of $n$ increases to $\mathcal{U}:=\{16,24\}$, the CPU time will increase to a considerable amount.

9 Conclusion

The core theoretical work of this paper locates in the uniform upper estimate for

[TABLE]

where $A^{(p)},A_{n}^{(p)}$ are defined in (2.4), (3.3) respectively. This determines the error estimate for approximation error and give a complete answer to regularization procedure.

In experiments, the algorithm has its advantage over other classical regularization method:

• •

It induces a noise-independent a-priori parameter choice strategy for function of a specific class

[TABLE]

where $p$ is the order of numerical differentiation. Good filtering effect (error approaches 0) is displayed when the algorithm acts on functions of above class with best parameter choice.

• •

Derivatives discontinuities can also be recovered well although there exists a unknown constant $\kappa$ to test in experiments.

Appendix A Representation of $M_{n}^{(p)}$

[TABLE]

where

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Appendix B Some Fourier expansions

Lemma B.1

For $A^{(p)}$ defined in (2.4), when $j\geq n+1$, set

[TABLE]

Then Fourier coefficients are determined as follows:

[TABLE]

[TABLE]

[TABLE]

Lemma B.2

For $A^{(p)}$ defined in (2.4), when $j\geq n+1$, set

[TABLE]

Then Fourier coefficients are determined as follows

[TABLE]

[TABLE]

[TABLE]

Appendix C Some Inequalities

Proposition C.1

Set $L_{n},K_{n},F_{n},T_{n}$ defined as in Theorem 4.1. Then

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Acknowledgement

I would like to thank Professor Nailin Du for his helpful discussions and suggestion that extend numerical differentiation into higher order. And the author thank Professor Hua Xiang for his comments. Besides, I really appreciate the help in figure modification from Nianci Wu. The work of Yidong Luo is supported by the National Natural Science Foundation of China under grants 11571265 and NSFC-RGC No. 11661161017.