A Fast Method for Numerical Realization of Fourier Tools

*Anry Nersessian*

### **Abstract**

This chapter presents new application of author's recent algorithms for fast summations of truncated Fourier series. A complete description of this method is given, and an algorithm for numerical implementation with a given accuracy for the Fourier transform is proposed.

**Keywords:** Fourier series, Fourier transforms, approximation, over-convergence phenomenon, numerical methods

### **1. Introduction**

One of the classical tools of Mathematics is the apparatus of Fourier series, based on the orthogonal system *<sup>e</sup><sup>i</sup>πkx* , *<sup>k</sup>* <sup>¼</sup> 0, � 1, � 2, … complete in *<sup>L</sup>*2½ � �1, 1 . However, in practice it is extremely limited because of the poor approximation of piecewise smooth functions. Thus, in the case when a function has discontinuity points (taking into account also discontinuities at the ends of the interval ½ � �1, 1 ), an intense oscillation arises in their neighborhood, and the uniform convergence is absent (the Gibbs phenomenon). This leads to a slow *L*2-convergence on the entire segment ½ � �1, 1 . The corresponding phenomenon is observed in the case of Fourier interpolation. Classical methods of summation truncated Fourier series (see, for example, [1], Chapter 3) do not actually change the situation.

Since the beginning of the last century, a huge amount of research has been carried out, devoted to the methods of effective summation of the Fourier series. Let us briefly dwell on the works to overcome this situation (for detiles see, for example, articles [2–13] and their links).

The pioneer of "overcoming the Gibbs phenomenon" is A.N.Krylov, who at the dawn of the 20th century (see [2]) proposed methods that he later developed in the monograph [3]. In particular, he proposed the following approach. Let a piecewise smooth function *f* is given on the segment ½ � �1, 1 with Fourier coefficients *fs* , *<sup>s</sup>* <sup>¼</sup> 0, � 1, … , � *n*, *n* ≥1, and with the following jumps of the function *f* or its derivatives until degrees of the order *q* ≥1 at the points f g *ak* , �1< *a*<sup>1</sup> < ⋯ < *am* ¼ 1, 1≤ *m* < ∞,

$$\begin{aligned} A\_{p,k}(f) &= f^{(k)}\left(a\_p - \mathbf{0}\right) - f^{(k)}\left(a\_p + \mathbf{0}\right), \\ k &= \mathbf{0}, \mathbf{1}, \dots, q \ge \mathbf{0}, p = \mathbf{1}, \dots, m. \end{aligned} \tag{1}$$

In the neighborhoods of other points we assume that *f* ∈*C<sup>q</sup>*þ<sup>1</sup> . Let us construct a function *g* ¼ *g x*ð Þ, *x*∈ ½ � �1, 1 with Fourier coefficients *gs* , 0 ≤ ∣*s*∣ ≤*n*, which has the same jumps at the same points, and *g* ∈*C<sup>q</sup>*þ<sup>1</sup> at the neighborhoods of other points.

Known the jumps (1), one can construct, e.g. piecewise-polynomial *g*. As a result, 2-periodic extension of the function *F* ¼ ð Þ *f* � *g* is *q* times continuously differentiable on whole axis, and then

$$f(\mathbf{x}) - \mathbf{g}(\mathbf{x}) = \sum\_{s=-n}^{n} (f\_s - \mathbf{g}\_s) e^{i\pi \mathbf{x} \mathbf{x}} + r\_s(\mathbf{x}),$$

where *rs*ð Þ¼ *<sup>x</sup> o s*�*<sup>q</sup>* ð Þ, *<sup>s</sup>* ! <sup>∞</sup>, *<sup>x</sup>*<sup>∈</sup> ½ � �1, 1 .

Therefore, taking into account only first 2ð Þ *n* þ 1 Fourier coefficients and truncating the remainder term *rs*, it is possible to approximate *f* in the form

$$f(\mathbf{x}) \simeq f\_n(\mathbf{x}) \stackrel{def}{=} \mathbf{g}(\mathbf{x}) + \sum\_{k=-n}^{n} (f\_s - \mathbf{g}\_s) e^{i\mathbf{x}\cdot\mathbf{x}} \tag{2}$$

to three decades, in the KE-method scheme, not only polynomials have been used as

The classical definition of the partial sums of the Fourier series is based on a gradual increase in the frequencies of the Fourier system. Below we use a more

are Fourier coefficients of *f*, and *Dn* ¼ f g *dk* , *k* ¼ 1, … , *n*, is a set of *n* different

*<sup>k</sup>P<sup>β</sup><sup>k</sup>* ð Þ *x* exp ð Þ *iπ λ<sup>k</sup> x* , where polynomials *P<sup>β</sup><sup>k</sup>* ð Þ *x* ≢ 0 have the exact degree *β<sup>k</sup>* and

*<sup>k</sup>* 1 þ *β<sup>k</sup>* ð Þ≤ *n*. The number *m* will be considered below as the degree of the

**Definition 3.** *We call the system* f g *λ<sup>k</sup>* ⊂ *as parameters of the quasi-polynomial q x*ð Þ¼ *<sup>k</sup>P<sup>β</sup><sup>k</sup>* ð Þ *x* exp ð Þ *iπ λ<sup>k</sup> x* ∈ *Qn and the number* 1 þ *β<sup>k</sup> as multiplicity of the parameter λk.* **Definition 4.** *Let parameters* f g *<sup>λ</sup><sup>k</sup>* , *<sup>k</sup>* <sup>¼</sup> 1, … , *n are* <sup>1</sup> *fixed. We denote by Qn*ð Þ f g *<sup>λ</sup><sup>k</sup>*

**Remark 1.** *The set of quasi-polynomials q x*ð Þ∈ *Qn is invariant with respect to a linear change of the variable x. The set of quasi-polynomials q x*ð Þ∈ *Qn*ð Þ f g *λ<sup>k</sup> is invariant*

For a fixed integer *n* ≥1, consider a set of non-integer parameters

<sup>1</sup> Unless otherwise stated, we will not exclude the repetition of values of *<sup>λ</sup><sup>k</sup>* in a set f g *<sup>λ</sup><sup>k</sup>* .

**Definition 2.** *Let n*≥1 *be a fixed integer. Consider a system of functions Un* ¼ f g exp ð Þ *iπ λ<sup>k</sup> x* , *λ<sup>k</sup>* ∈ , *x*∈½ � �1, 1 , *k* ¼ 1, 2, … , *n, where* f g *λ<sup>k</sup> are arbitrary parameters.*

*Consider the linear span Qn* ¼ *span U*f g*<sup>n</sup> . We call a function q*∈ *Qn as quasi-*

*f <sup>k</sup>* exp ð Þ *iπ kx* , *x*∈½ � �1, 1 , (7)

�*<sup>i</sup> <sup>π</sup> s x dx*, *<sup>s</sup>* <sup>¼</sup> 0, � 1, � 2, … (8)

, if and only if when either *q x*ð Þ� 0 or *q x*ð Þ¼

**Definition 1.** *Let us call the truncated Fourier series any sum of the form*

with a complete description of the required algorithms from [17].

**2. Acceleration of convergence of Fourier series**

*A Fast Method for Numerical Realization of Fourier Tools*

*DOI: http://dx.doi.org/10.5772/intechopen.94186*

*Sn*ð Þ¼ *x*

*fs* <sup>¼</sup> <sup>1</sup> 2 ð1 �1 *f x*ð Þ*e*

integers (*n*≥1). We will assume that *D*<sup>0</sup> ¼ Ø.

*polynomial of degree at most n.* It is easy to see that *<sup>q</sup>*<sup>∈</sup> *Qn* <sup>P</sup>

*the set of all corresponding q x*ð Þ∈ *Qn.*

*with respect to a shift of the variable x.*

f g *λ<sup>k</sup>* ⊂ , *k*∈ *Dn*, and the following infinite sequence

*def* X *k* ∈ *Dn*

In [14], for a smooth function *f*, a combination of the solution of Eq.(6) and the Pade approximation (applied to the asymptotic expansion of the Fourier coefficients) is implemented there. As a result of numerical experiments it turned out that, as a rule, such an algorithm is "almost exact" on certain infinite-dimensional spaces, although only a finite number of Fourier coefficients is used. In [15], a similar phenomenon was discovered for Bessel series, and in [16] for some eigenfunction expansions. The theoretical substantiation of this phenomenon for Fourier series and the corresponding algorithm were published in [17–19]. Our goal is to use this method in the case of one-dimensional Fourier transforms. Let us start

functions of *g*.

**2.1 Basic definitions**

where

*<sup>m</sup>* <sup>¼</sup> <sup>P</sup>

P

**31**

quasi-polynomial *q*.

general notations from the work [17].

However, the computing of the jumps f g *Ask*ð Þ*f* directly by the function *f* sharply limits the scope of practical application of the method of A. N. Krylov.

The "spectral" method proposed by Kurt Eckhoff in [4] (1993) turned out to be more practical, as it is based only on the use of coefficients *fs* � � (see also [5]). It is easy to obtain using the integration in parts the following asymptotic representation of the Fourier coefficients:

$$f\_s = -\frac{1}{2} \sum\_{p=1}^{m} e^{-i\pi s d\_p} \sum\_{k=0}^{q-1} \frac{A\_{p,k}}{(i\pi s)^{k+1}} + r\_s,\\ r\_n = o(s^{-q}), s \to \infty. \tag{3}$$

As a function *g* from (2) K. Eckhoff used those Bernoulli polynomials f g *Bk*ð Þ *x* , *x*∈ ½ � �1, 1 , *k*≥0, which Fourier coefficients f g *bk*,*<sup>s</sup>* have the following simple form

$$b\_{k, \mathfrak{s}} = \begin{cases} 0, & \mathfrak{s} = \mathbf{0}, \\ \frac{(-1)^{s+1}}{2\left(i\pi\mathfrak{s}\right)^{k+1}}, & \mathfrak{s} \neq \mathbf{0}, k = \mathbf{1}, 2, \dots \end{cases} \tag{4}$$

Denoting *B*0ð Þ¼ *x* 1, polynomials f g *Bk*ð Þ *x* , *k* ¼ 0, 1, … , *n*, *x*∈½ � �1, 1 compose a basis on the space of polynomials of degree *n*. Bernoulli polynomials extended to the real axis with period 2 as piecewise-smooth functions.

According to Krylov's scheme, the sequence

$$F\_n(\mathbf{x}) \stackrel{def}{=} \sum\_{p=1}^m \sum\_{k=0}^q A\_{p,k} B\_k(\mathbf{x} - a\_p - \mathbf{1}) + \sum\_{\varepsilon=-n}^n a\_\varepsilon e^{i\mathbf{x}\cdot\mathbf{x}},\tag{5}$$

where the quantities f g *ω<sup>s</sup>* are given explicitly, converges to *f* with the rate *o n*�*<sup>q</sup>* ð Þ, *<sup>n</sup>* ! <sup>∞</sup>.

K. Eckhoff suggests to find approximate values of jumps *A*~ *<sup>p</sup>*,*<sup>k</sup>*≈*Ap*,*<sup>k</sup>* � � by solving the following system of linear equations with the Vandermonde matrix, obtained by principal part (1) choosing the indexes *s* ¼ *sk*, *k* ¼ 1, 2, … , *m q*ð Þ þ 1 , *θn*≤ ∣*sk*∣ ≤*n*, 0<*θ* ¼ *const* <1

$$f\_s = -\frac{1}{2} \sum\_{p=1}^{m} e^{-i\pi s d\_p} \sum\_{k=0}^{q} \frac{\tilde{A}\_{p,k}}{\left(i\pi s\right)^{k+1}}, s = s\_1, s\_2, \dots, s\_{m(q+1)}\tag{6}$$

Further, the problem is solved by the Krylov method. It is natural to call this acceleration scheme the Krylov-Eckhoff method **(KE-method)**. Over the past two *A Fast Method for Numerical Realization of Fourier Tools DOI: http://dx.doi.org/10.5772/intechopen.94186*

to three decades, in the KE-method scheme, not only polynomials have been used as functions of *g*.

In [14], for a smooth function *f*, a combination of the solution of Eq.(6) and the Pade approximation (applied to the asymptotic expansion of the Fourier coefficients) is implemented there. As a result of numerical experiments it turned out that, as a rule, such an algorithm is "almost exact" on certain infinite-dimensional spaces, although only a finite number of Fourier coefficients is used. In [15], a similar phenomenon was discovered for Bessel series, and in [16] for some eigenfunction expansions. The theoretical substantiation of this phenomenon for Fourier series and the corresponding algorithm were published in [17–19]. Our goal is to use this method in the case of one-dimensional Fourier transforms. Let us start with a complete description of the required algorithms from [17].
