Abstract
In this article we consider the numerical differentiation of periodic functions specified by noisy data. A new method, which is based on the truncated singular value decomposition (TSVD) regularization technique of a suitable compact operator, is presented and analysed. It turns out the new method coincides with some type of truncated Fourier series approach. Numerical examples are also given to show the efficiency of the method.
1. Introduction
Numerical differentiation is the problem of determining derivatives of a function from its values on an interval or some scattered points. It arises from many scientific researches and applications, e.g. the identification of the discontinuous points in an image process Citation1; the problem of solving the Abel integral equation Citation2; the problem of determining the peaks in chemical spectroscopy Citation3, some inverse problems in mathematical physics equations Citation4, etc. The main difficulty is that it is an ill-posed problem, which means small errors in the measurements of a function may lead to large errors in its computed derivatives Citation4,Citation5. A number of techniques have been developed for numerical differentiation Citation4–11. In Citation10, a theoretical result has been obtained for arbitrary numerical differentiation of a function in Hn[0, 1]. But in their method, different solution processes are used for different n. In this article, we will give a method for numerical differentiation of periodic functions. In our method, the solution processes will be uniform for different n, which means that the method is self-adaptive.
One type of method of the numerical differentiation problem is to transform it into an operator equation of the first kind. In fact, for given g ∈ H1[0, 1], finding f = Dg = g′ is equivalent to solving the Volterra integral equation of the first kind
(1)
In Citation11, a modified operator of K1 in (1) has been presented for general functions. In this article, another modified operator will be considered for periodic functions. For the new operator, better properties will be obtained and we can successfully deal with high-order differentiation, as well as first order. In the present situation a singular system of the new operator can be obtained easily. Therefore, it seems natural to use the truncated singular value decomposition (TSVD) method for our problem.
The article is organized as follows. In Section 2, we will discuss a singular system and the properties of the new operator. In Section 3, an algorithm based on the TSVD method for first-order numerical differentiation is proposed and analysed, and in Section 4 it is extended to high-order differentiation. Some numerical results are given in Section 5, which show that the new method works well.
2. Singular system and its properties
Some notations will be used in the sequel. Let ℂ1 = ℂ1(ℝ) denote the space of continuous functions with period 1, the associated norm is given by
(2)
and let
denote the space of square integrable functions on the interval [0, 1] with period 1, the associated norm is given by
(3)
Finally, let
It is obvious that the operator K1 is the inverse of the following operator T1:
(4)
We consider now a slightly different operator T defined by
(5)
where
In this article, ⟨·, ·⟩ denotes the inner product in space L2. Then T is invertible and K = T−1 : is given by
(6)
K is compact and its adjoint operator K* = −K. It can be verified that the eigenproblem K*Ke = λe is equivalent to
(7)
Thus the eigenvalues of K*K are
There are exactly two eigenfunctions corresponding to each eigenvalue λl:
Therefore the singular values and the corresponding singular functions of operator K can be taken as
(8)
For a function
, the solution to the equation
(9)
can be expressed by the singular system of the operator K as follows Citation5:
(10)
The above formula shows that the expansion of function f by the singular functions
is just the Fourier series of f on [0, 1]. The important thing is that these Fourier coefficients
can be calculated in terms of the singular values
and the Fourier coefficients
of the given function g.
In the following, we discuss the generalized smooth scale of a function with respect to the operator K, which is a key index concerning convergence rates of a regularization method Citation5. We need some results of Fourier series. For a periodic function φ(t) with period 2l, its Fourier series has the form
or
where
LEMMA 2.1 Citation12
Suppose that function φ(t) is 2l-periodic and of bounded variation on [0, 2l], and let be the Fourier coefficients of φ(t). Then
(11)
Suppose now , then from integration by parts we have
(12)
and
(13)
On the other hand, we know that for a function
,
, if and only if Citation4
(14)
Thus from Lemma 2.1, we have the following theorem.
THEOREM 2.2
Suppose , then we have
(15)
Moreover, if Dn f is of bounded variation, then
(16)
Especially if
, then
(17)
Proof
If , let
. If n is even, then we have
(18)
If n is odd, then we let
. Then h1 = −K*h, and therefore
(19)
It is well known that
(20)
Hence (14) holds. Moreover, if Dn f is of bounded variation, then (15) results from Lemma 2.1 and (12–13).▪
Remark 1
Generally, a function is not in
. Let
(21)
Then
and
. Thus the differentiation problem for g can be transferred to solving the following operator equation
(22)
Remark 2
For a non-periodic function g ∈ Hn[0, 1](e.g. g(t) = exp(t)), let be the periodic continuation of function g from the interval [0, 1] to all ℝ with period 1. Generally, we have
, even if we let
(23)
we only have
. Therefore, the present method can only be used in the following when the periodic continuation does not break the smoothness of the function.
3. First-order differentiation
Instead of g(s), in practice we usually only have perturbed data gδ(s). Furthermore, the condition
(24)
is assumed with δ being a known error level. The perturbed form of Equation (21) is
(25)
where
(26)
We have
(27)
and it is clear that
. The TSVD solution
to Equation (24) can be given as follows:
(28)
where m = m(δ) is determined by the discrepancy principle
(29)
From the theory of TSVD method we have the following theorem.
THEOREM 3.1 Citation13
Let gδ satisfy (23), τ > 1, is defined by (27) and (28). Suppose
, ν > 0, then
(30)
and
(31)
where
(32)
From Theorem 3.1 and Remark 1, we have the following result.
COROLLARY 3.2
Suppose ‖gδ − g‖ ≤ δ, f = g′, , then we have
(33)
Moreover, if Dn g is of bounded variation, then for
, the TSVD solution
defined above satisfies
(34)
4. High-order differentiation
Let
(35)
and the operator
is compact. For a function
, its k-th derivative h = Dkg(1 ≤ k ≤ n) satisfies the equation
(36)
The perturbed version of Equation (35) is
(37)
Similar to the discussion in the above, the singular values and the corresponding singular functions of operator Bk can be taken as
(38)
If k is even,
(39)
If k is odd,
(40)
Analogously, we have the following theorem.
THEOREM 4.1
Suppose , then
(41)
Moreover, if Dn h is of bounded variation on [0, 1], then
(42)
Especially, if
, then
(43)
The TSVD solution to Equation (36) is given by
(44)
where m = m(δ) is determined by the discrepancy principle
(45)
Remark 3
We have
(46)
Therefore, for 1 < k ≤ n, the value of m in (43) is equal to the one in (28).
We have the following result.
THEOREM 4.2
Let gδ satisfies (23), τ > 1, is defined by (42) and (43), Suppose
, ν > 0, then
(47)
and
(48)
Especially, we have the following corollary.
COROLLARY 4.3
Suppose , h = Dkg, gδ satisfies (23), τ > 1,
is defined by (42) and (43), then we have
(49)
Moreover, if h is of bounded variation on [0, 1], then
the TSVD solution hδ defined above satisfies
(50)
Next, we will give a uniform convergence result.
LEMMA 4.4 Citation14
If h ∈ H1[0, 1], then
(51)
where C0 is a constant independent of f.
It can be verified that if is the TSVD solution to the equation
, then
is the TSVD solution to the equation
. Thus from Corollary 4.3 and Lemma 4.4, we have the following theorem.
THEOREM 4.5
Suppose and Dng is of bounded variation on [0, 1], and the estimate (23) is fulfilled, then
converges uniformly to Dkg as δ → 0, and
,
(52)
5. Numerical examples
In this section, we present numerical results of some test examples. In all cases, the discretization knots ti = ih, i = 1, …, N, with N = 2048, h = 1/N. The perturbed discrete data is given by
where
are generated by function (2 × rand(N + 1, 1) − 1) × δ1 in Matlab.
In computation, one has to calculate the integrals
First we obtain the cubic spline interpolation function s(t) for non-oscillatory function gδ(t) with the periodic boundary conditions
In this case, the function s(t) can be regarded as the perturbed data of function g(t) and the integrals
where Mi = s″(ti). We only need the error estimates δ in
As in our example the error satisfies a uniform distribution, we let
in practical computing.
All examples are computed by using Matlab with the parameter τ = 1.01. The L2-norm relative errors and the maximum-norm relative errors are given in and to verify the theoretical results. and with δ1 = 0.5 and 0.1 compare qualitatively the computed solutions and the exact ones.
Let us first consider the numerical comparison between the present method (M1) and the method in Citation11 (M2) for computing the first-order derivative.
Here and in what follows, the solid curve means the exact solution, while the dotted–dashed curve indicate the solutions of M1 and M2, respectively.
Example 1
First we consider the periodic function
From Corollary 3.2 and Theorem 3.1 in Citation11, we have that the convergence rates of M1 and M2 are O(δ) and O(δ3/5), respectively. The numerical results are shown in and .
The numbers in columns 3, 4 and 7, 8 show that the results of M1 are much better than M2 for periodic function. Moreover, the numbers in columns 5 and 9 exhibit that the ratios are stable when δ → 0.
Table 1. Comparison of M1 and M2 for periodic function.
Example 2
We consider the non-periodic function
From Remark 2, Corollary 3.2 and Theorem 3.1 in Citation11, we have that the convergence rates of M1 and M2 are O(δ1/3) and O(δ3/5), respectively. The numerical results are shown in and . The numbers in columns 3, 4 and 7, 8 show that the results of M1 are inferior to the ones of M2.
These two examples show that the present method is fit for solving the derivatives of periodic smooth functions.
Next, we give the numerical results of the high-order derivatives of Example 1 in and . All results show that the new method works well.
Figure 3. Results of high-order derivatives of Example 1: (a) second-order derivative with δ1 = 0.1 and (b) third-order derivative with δ1 = 0.1.
![Figure 3. Results of high-order derivatives of Example 1: (a) second-order derivative with δ1 = 0.1 and (b) third-order derivative with δ1 = 0.1.](/cms/asset/48937e37-5c31-4f1e-9b5f-3dd93cf7c835/gipe_a_492517_f0003.gif)
Table 2. Comparison of M1 and M2 for non-periodic function.
Table 3. Results of high-order derivatives of Example 1.
Example 3
Citation8 We consider the function
Let
be the periodic continuation of g from the interval [0, 1] to all ℝ with period 1. Thus the original problem can be changed to find derivatives of
. Since g(k)(0) = g(k)(1) = 0, ∀k ∈ ℕ, then
. From the Corollaries 3.2 and 4.3, we have that the convergence rate of the present method is O(δ). The numerical results are shown in and and .
We conclude the article by giving some remarks.
1. | A new method for the numerical differentiation of periodic functions is proposed and analysed. All the test numerical examples presented in the article show that the new method works well. | ||||
2. | The new method proposed in the article can be regarded as a truncated Fourier series method. The important thing is that the new operators K build a bridge between Fourier truncation method and regularization method. So we can choose a regularization parameter by a posteriori rule. |
Figure 4. Results of first-order derivative of Example 3: (a) first-order derivative with δ1 = 0.5 and (b) first-order derivative with δ1 = 0.1.
![Figure 4. Results of first-order derivative of Example 3: (a) first-order derivative with δ1 = 0.5 and (b) first-order derivative with δ1 = 0.1.](/cms/asset/b86851f3-4a3f-4ec7-8d91-25a03c1c82b4/gipe_a_492517_f0004.gif)
Figure 5. Results of high-order derivatives of Example 3: (a) second-order derivative with δ1 = 0.1 and (b) third-order derivative with δ1 = 0.1.
![Figure 5. Results of high-order derivatives of Example 3: (a) second-order derivative with δ1 = 0.1 and (b) third-order derivative with δ1 = 0.1.](/cms/asset/4b025ca2-f7f7-42e8-ab67-5290b31f545f/gipe_a_492517_f0005.gif)
Table 4. Results of the first-order derivative of Example 3.
Acknowledgements
The authors are grateful to Prof. Daniel Lesnic and anonymous referees for their valuable suggestions that improved the article substantially. The project is supported by the Science Foundation of Guangdong Ocean University, No. 0812279.
References
- Deans, SR, 1983. Radon Transform and its Applications. New York: John Wiley & Sons; 1983.
- Gorenflo, R, and Vessella, S, 1991. Analysis and Applications of Abel Integral Equations, Lecture Notes in Mathematics. Vol. 1461. Berlin: Springer; 1991.
- Hegland, M, 1998. Resoluton Enhancement of Spectra Using Differentiation. Australian National University; 1998, preprint.
- Hanke, M, and Scherzer, O, 2001. Inverse problems light: Numerical differentiation, Amer. Math. Monthly 108 (2001), pp. 512–521.
- Engl, HW, Hanke, M, and Neubauer, A, 1996. Regularization of Inverse problems. Dordrecht: Kluwer Academic; 1996.
- Murio, DA, 1993. The Mollification Method and Numerical Solution of Ill-posed Problems. New York: John Wiely & Sons; 1993.
- Murio, DA, Mejia, CE, and Zhan, S, 1998. Discrete mollification and automatic numerical differentiation, Comput. Math. Appl. 35 (1998), pp. 1–13.
- Qian, Z, Fu, CL, Xiong, XT, and Wei, T, 2006. Fourier truncation method for high order numerical derivatives, Appl. Math. Comput. 181 (2006), pp. 940–948.
- Ramm, AG, and Smirnova, AB, 2001. On stable numerical differentiation, Math. Comput. 70 (2001), pp. 1131–1153.
- Wang, YB, Hon, YC, and Cheng, J, 2006. Reconstruction of high order derivatives from input data, J. Inverse Ill-posed Probl. 14 (2006), pp. 205–218.
- Zhao, ZY, Meng, ZH, and He, GQ, 2009. A new approach to numerical differentiation, J. Comput. Appl. Math. 232 (2009), pp. 227–239.
- Zygmund, A, 2002. Trigonometric Series. Cambridge, UK: Cambridge University Press; 2002.
- Vainikko, GM, 1982. The discrepancy principle for a class of regularization methods, USSR Comput. Math. Phys. 22 (1982), pp. 1–19.
- Adams, RA, 1975. Sobolev Spaces, Pure and Applied Mathematics. Vol. 65. New York: Academic Press; 1975.