Anry Nersessian
Institute of Mathematics, National Academy of Sciences of the Republic of Armenia, Yerevan, Armenia.
DOI: 10.4236/apm.2024.147033   PDF    HTML   XML   142 Downloads   732 Views   Citations

Abstract

In the author’s recent publications, a parametric system biorthogonal to the corresponding segment of the exponential Fourier system was unusually effective. On its basis, it was discovered that knowledge of a finite number of Fourier coefficients of function f from an infinite-dimensional set of elementary functions allows f to be accurately restored (the phenomenon of over-convergence). Below, parametric biorthogonal systems are constructed for classical trigonometric Fourier series, and the corresponding phenomena of over-convergence are discovered. The decisive role here was played by representing the space L₂ as an orthogonal sum of two corresponding subspaces. As a result, fast parallel algorithms for reconstructing a function from its truncated trigonometric Fourier series are proposed. The presented numerical experiments confirm the high efficiency of these convergence accelerations for smooth functions. In conclusion, the main results of the work are summarized, and some prospects for the development and generalization of the proposed approaches are discussed.

Keywords

Fourier Series, Acceleration of Convergence, Parametric Biorthogonalization, Spectral Methods, Over-Convergence Phenomenon

Share and Cite:

1. Introduction

Let’s start with an excursion into the two-century history of studying the classical Fourier series, stopping only at some of its stages. We do not claim to be complete about the review as a whole or the exhaustive nature of the selected stages. The scientists mentioned below (perhaps subjectively) are just a few who have made significant contributions to the topic of this work, which can be described as the problem of approximate restoration of a function defined on a finite segment using only a finite set of its Fourier coefficients.

One can read about our former fast algorithms based on the application of the phenomenon of super-convergence for the exponential Fourier series and the goals of this work in Sections 1.5 and 1.6.

1.1. The Brilliant Legacy of Joseph Fourier

One of the classical tools of mathematics is the apparatus of Fourier series based on the orthogonal system $\left\{{\text{e}}^{i\pi kx}\right\},k=0,\pm 1,\pm 2,\cdots$ , complete in $L_2\left[-\mathrm{1,1}\right]$ . The straight summation method for truncated Fourier series is the following

$f\left(x\right)\simeq {\displaystyle \sum _{k=-n}^n}{f}_k{\text{e}}^{i\pi kx}\mathrm{,}{f}_k=\frac{1}{2}{\displaystyle {\int }_{-1}^1}f\left(t\right){\text{e}}^{-i\pi kx}\text{d}t\mathrm{,}x\in \left[-\mathrm{1,1}\right]\mathrm{,}$ (1.1)

where $n\ge 0$ is an integer and $f\left(x\right)\in L_2\left[-\mathrm{1,1}\right]$ .

Expressing the analogical system in terms of sines and cosines, Jean-Baptiste Joseph Fourier (1768-1830) first used it in the early 19th century to approximately solve the heat equation. The problem he posed of representing a general function in a series according to this system caused one of the most significant scientific discussions in natural science and formed the basis for the further development of physics and mathematics.

Fourier series and their numerous generalizations played an essential role in solving many theoretical and applied problems until the mid-20th century. However with the advent of computers and the use of such new means, such as finite difference methods and wavelets, the practical possibilities of the Fourier series have faded into obscurity. It turned out that when approximating a smooth function $f\left(x\right)$ by Formula (1.1), in the vicinity of the ends of segment [−1, 1], intense oscillations arise, and there is no uniform convergence at $n\to \infty$ (the Gibbs phenomenon¹). As a rule, even the L₂-convergence of approximation (1.1) to a smooth function $f\left(x\right)$ $n\to \infty$ is very slow. Other classical methods of summation of Fourier series only slightly speeded the convergence for smooth functions (see, for example, [1], Volume 1, Chapter III). For a long time, it seemed that the Fourier series’ applied capabilities had been exhausted.

1.2. The Approach of Alexey N. Krylov

The pioneer of “overcoming the Gibbs phenomenon for Fourier series” was undoubtedly an outstanding scientist and shipbuilding engineer A. N. Krylov, who, even at the dawn of the twentieth century (see [2], 1907), proposed numerical methods, which later was included in his classical book [3] (1911). In particular, he proposed the following approach.

Let a piece-wise smooth function f is given on the segment $\left[-\mathrm{1,1}\right]$ with Fourier coefficients $\left\{f_s\right\}$ , $s=\mathrm{0,}\pm \mathrm{1,}\cdots \mathrm{,}\pm n$ , $n\ge 1$ , and with the following jumps of the function f and its derivatives until degrees of the order $q\ge 1$ at the points $\left\{a_k\right\}$ , $-1=a_1<\cdots <a_m=1$ , $1\le m<\infty$ ,²

$A_{p\mathrm{,}k}\left(f\right)=f^{\left(k\right)}\left(a_p-0\right)-f^{\left(k\right)}\left(a_p+0\right)$ (1.2)

where $q\ge 0$ , $k=\mathrm{0,1,}\cdots \mathrm{,}q$ , $p=\mathrm{1,}\cdots \mathrm{,}m$ .

In the neighborhoods of other points, we assume that $f\in C^{q+1}$ . Let us construct a function $g=g\left(x\right)$ , $x\in \left[-\mathrm{1,1}\right]$ , with Fourier coefficients $\left\{g_s\right\}$ , $0\le \left|s\right|\le n$ , which has the same jumps at the same points, and $g\in C^{q+1}$ at the neighborhoods of other points. Known as the jumps (1.2), one can construct, e.g., piece-wise polynomial g. As a result, a 2-periodic extension of the function $F=\left(f-g\right)$ is q times continuously differentiable on the whole x-axis. Therefore, taking into account only the first $\left(2n+1\right)$ Fourier coefficients and truncating the remainder term $r_s$ , it is possible to approximate f in the form

$f\left(x\right)\simeq F_n\left(x\right)=g\left(x\right)+{\displaystyle \sum _{k=-n}^n}\left(f_s-g_s\right){\text{e}}^{i\pi sx}\mathrm{.}$ (1.3)

Unfortunately, the ideas of A. N. Krylov remained in the shadow for a long time in the international arena, and such a method became widely known thanks mainly to the works of Cornelius Lanczos, which were published half a century later (see [4], 1966).

1.3. Spectral Method of Knut S. Eckhoff

However, the computing of the jumps $\left\{A_{sk}\left(f\right)\right\}$ directly by the function f extremely limits the scope of the practical application of Krylov’s method (1.3).

K. S. Eckhoff developed in [5] (1993) and [6] the “spectral” method, which turned out to be much more practical since it is based only on the use of Fourier coefficients $\left\{f_s\right\}$ . Let us briefly describe this approach.

By using the integration in part, it is easy to obtain the following asymptotic representation of the Fourier coefficients:

$f_s=-\frac{1}{2}{\displaystyle \sum _{p=1}^m}{\text{e}}^{-i\pi s{a}_p}{\displaystyle \sum _{k=0}^{q-1}}\frac{A_{p,k}}{{\left(i\pi s\right)}^{k+1}}+r_s,r_n=o\left(s^{-q}\right),s\to \infty .$ (1.4)

As a function of g from (1.3) K. Eckhof used Bernoulli polynomials $\left\{B_K\left(x\right)\right\}$ , $x\in \left[-\mathrm{1,1}\right]$ , $k\ge 0$ . Their Fourier coefficients $\left\{b_{k\mathrm{,}s}\right\}$ have the following simple form

$b_{k,s}=(\begin{array}{ll}0,\hfill & s=0,\hfill \\ \frac{{\left(-1\right)}^{s+1}}{2{\left(i\pi s\right)}^{k+1}},\hfill & s\ne 0,k=1,2,\cdots .\hfill \end{array}$

Denoting $B_0\left(x\right)=1$ , polynomials $\left\{B_k\left(x\right)\right\}$ , $k=\mathrm{0,1,}\cdots \mathrm{,}n$ , $x\in \left[-\mathrm{1,1}\right]$ , are composed of a basis on the space of polynomials of degree n. Extended to the real axis with period 2, Bernoulli polynomials are piece-wise smooth functions. This polynomial can be calculated recursively as following

$B_0\left(x\right)=\mathrm{1,}{B}_k\left(x\right)={\displaystyle \int B_{k-1}\left(x\right)\text{d}x}\mathrm{,}{\displaystyle {\int }_{-1}^1}{B}_k\left(x\right)\text{d}x=\mathrm{0,}k=\mathrm{2,3,}\cdots \mathrm{.}$

According to Krylov’s scheme (1.3),

$f\left(x\right)\simeq F_n\left(x\right)={\displaystyle \sum _{p=1}^m}{\displaystyle \sum _{k=0}^q}{A}_{p,k}{B}_k\left(x-a_p-1\right)+{\displaystyle \sum _{s=-n}^n}{\omega }_s{\text{e}}^{i\pi sx}\mathrm{,}$ (1.5)

where the quantities $\left\{{\omega }_s\right\}$ are given explicitly, converges to f with the rate $o\left(n^{-q}\right)\mathrm{,}n\to \infty$ .

K. Eckhoff suggests finding approximate values of jumps $\left\{{\tilde{A}}_{p\mathrm{,}k}\approx A_{p\mathrm{,}k}\right\}$ by solving the following system of linear equations with the Vandermonde matrix, which is obtained by the principal part of (1.4) choosing the indexes $s=s_k$ , $k=1,2,\cdots ,m\left(q+1\right)$ , $\theta n\le \left|s_k\right|\le n$ , $0<\theta =const<1$ .

$f_s=-\frac{1}{2}{\displaystyle \sum _{p=1}^m}{\text{e}}^{-i\pi s{a}_p}{\displaystyle \sum _{k=0}^q}\frac{{\tilde{A}}_{p,k}}{{\left(i\pi s\right)}^{k+1}},s=s_1,s_2,\cdots ,s_{m\left(q+1\right)}.$(1.6)

We call this algorithm of acceleration Krylov-Eckhoff method (KE-method).

1.4. Briefly about Some Development of This Topic

The selection of works [7]-[13], together with their literature, to a certain extent, reflects “overcoming the Gibbs phenomenon” after A. Krylov’s research (see Section 1.2 above) and before K. Eckhoff’s spectral start. Here, first of all, it is necessary to emphasize the role of Cornelius Lanczos, whose works cover both theoretical and applied aspects of the Fourier series. His classic work has recently been reissued and, through the efforts of J. Boyd, significantly expanded (see [7]). In their work [8], James Geer and Nana S. Banerjee applied A. Krylov’s method, choosing an original system of trigonometric functions as the function g (see (1.3)). In [9] David Gottlieb and Chi-Wang Shu use the Fourier coefficients multiplied by a filter $\sigma \left(k/n\right)$ before calculating the jumps of the function $f\left(x\right)$ (see also [9] [10]). Article [11] of Herbert H. H. Homeier studies and generalizes a Levin-type algorithm for accelerating the Fourier series’s convergence involving a frequency parameter. In his work [12], some methods that work in the vicinity of jumps and other singularities of the Fourier series are presented.

The intensity of research has increased significantly in the 21st century, largely due to the use of the KE-method and its modifications. When studying Fourier interpolation based on a truncated FFT, it is shown in [13] that the Gibbs phenomenon does not disappear, but its “magnitude” can be effectively reduced by several decimal orders. Tobin A. Driscoll and Bengt Fornberg used in [14] a Padé-based algorithm for overcoming the Gibbs phenomenon.

In the article [15] of Dmitry Batenkov, a reconstruction algorithm for piecewise-smooth functions with a priori known smoothness and a number of jumps from their Fourier coefficients is provided, possessing the maximal possible asymptotic rate of convergence (see also [16]). The articles [17]-[20] of Ben Adcock are devoted to various aspects of “overcoming the Gibbs phenomenon”. In the work [21] of Yun Beong In, Krylov’s method with a weight function is used. In article [22] of Artur Barkhudaryan, Rafael Barkhudaryan, and Arnak Poghosyan, the asymptotic behavior of a function’s jumps is studied, and the exact asymptotic L₂ constants of the rate of convergence of the KE-method are computed. In [23], a nonlinear method for accelerating the convergence of the Fourier series was developed based on the application of Padé approximants to the asymptotic formula for Fourier coefficients. A recent tutorial [24] provides a detailed review and discussion of the topic of the Fourier series.

1.5. Over-Convergence Phenomenon for Fourier Series

New approaches have been demonstrated in works [25]-[30]. The following scheme was proposed.

Definition 1. We call any sum of the form

$S_n\left(x\right)={\displaystyle \sum _{k\in D_n}}{f}_k\exp \left(i\pi kx\right)\mathrm{,}x\in \left[-\mathrm{1,1}\right]\mathrm{,}$ (1.7)

the truncated Fourier series, where $D_n=\left\{d_k\right\}$ , $k=\mathrm{1,}\cdots \mathrm{,}n$ , is a set of n different integers ($n\ge 2$ ).

Definition 2. Let $n\ge 1$ be a fixed integer. Consider a system of functions $U_n=\left\{\exp \left(i\pi {\lambda }_{k}x\right)\right\}$ , ${\lambda }_k\in ℂ$ , $x\in \left[-\mathrm{1,1}\right]$ , $k=\mathrm{1,2,}\cdots \mathrm{,}n$ , where $\left\{{\lambda }_k\right\}$ are arbitrary parameters. Consider the linear span $Q_n=span\left\{U_n\right\}$ . We call a function $q\in Q_n$ as a quasi-polynomial of degree at most n.

Note that for a $h\ne 0$ , $h^{-1}\left(\exp \left(i\pi \left({\lambda }_k+h\right)x\right)-\exp \left(i\pi {\lambda }_{k}x\right)\right)\to i\pi x\exp \left(i\pi {\lambda }_{k}x\right)$ as $h\to 0$ . Therefore, it is clear that $q\in Q_n$ if and only if either $q\left(x\right)\equiv 0$ or $q\left(x\right)={\displaystyle {\sum }_k}{P}_{{\beta }_k}\left(x\right)\exp \left(i\pi {\lambda }_{k}x\right)$ , where the polynomials $P_{{\beta }_k}\left(x\right)\cancel{\equiv }0$ are of degree exactly ${\beta }_k$ , and $m={\displaystyle {\sum }_k}\left(1+{\beta }_k\right)\le n$ . The number m will be understood below as the degree of the quasi-polynomial q.

To properly accelerate the convergence of the Fourier series, a new “adaptive” nonlinear algorithm $\mathfrak{A}$ was applied, which contains the pseudo-inversion of one $m\times m$ matrix and the calculation of all zeros of a m-degree polynomial (for some details, see Sections 3.6.1 below). The following result was proven.

Theorem 1 (The phenomenon of the over-convergence [25] [27]). Let $f\in Q_m$ , the sets $D_m$ , ${\tilde{D}}_m$ , $D_n\cap {\tilde{D}}_m=\varnothing$ and the Fourier coefficients $\left\{f_s\right\}$ , $s\in D_m\cup {\tilde{D}}_m$ , of the function f be given. Denote by Λ the set of integer parameters in the approximation $f\simeq {\displaystyle {\sum }_k}{P}_{m_k}\left(x\right)\exp \left(i\pi {\mu }_{k}x\right)$ . In order for the approximation by Algorithm $\mathfrak{A}$ to be exact (that is, $f\left(x\right)\equiv F_m\left(x\right)$ , $x\in \left[-\mathrm{1,1}\right]$ ), it is necessary and sufficient that $\Lambda \subseteq D_m\cup {\tilde{D}}_m$ .

Unexpectedly, it turned out that the partial sum of the Fourier series is theoretically exact in the infinite-dimensional space ${\tilde{D}}_m$ and not only on the 2 m-dimensional set of exponentials $\exp \left(i\pi sx\right)$ with integer values of s. A huge difference! Of course, this led to a super-fast approximation of a smooth function F using its truncated Fourier series. Even more unexpectedly, that a similar over-convergence phenomenon was discovered for Fourier interpolation (i.e., for “finite Fourier series,” see [28])

1.6. The Purpose of This Work

We study the generalization of the adaptive algorithm of work [29] to the cases of trigonometric Fourier series. Initially, we hoped to follow the recommendations of works [29] and [30]. However, the case of the exponential Fourier series turned out to be unique (see previous section). The phenomenon of over-convergence for sine and cosine Fourier series is discovered below only after the corresponding parallelization of the acceleration algorithms.

The theoretical justification of our algorithms is based on the analytical method and, due to the large number of formulas, is concentrated in Sections 2 and 3. Numerical results are given in Section 4. The conclusion (Section 5) contains some comments on the presented results and a short discussion of possible generalizations to the cases of expansions by eigenfunctions of selfadjoint boundary value problems for general even-order ODEs with smooth coefficients on a finite segment.

2. Preliminary Formulas

The formulas below sometimes contain division by zero. We will present them formally for now and eliminate the singularities later as we go.

The main result here is the construction of corresponding parametric biorthogonal systems (see Section 2.3).

2.1. Fourier Cosine Series

For $f\in L_2\left[\mathrm{0,1}\right]$ , the following expansion is considered

$f\left(x\right)={\displaystyle \sum _{k=0}^{\infty }}{f}_k\cos \left(\pi kx\right),x\in \left[0,1\right],$ (2.1)

where

$f_0={\displaystyle {\int }_0^1}f\left(t\right)\text{d}t,f_s=2{\displaystyle {\int }_0^1}f\left(t\right)\cos \left(\pi kt\right)\text{d}t,s=1,2,\cdots .$

Let’s denote for an integer $n\ge 0$

$\left(n_1\mathrm{,}{n}_2\right)=(\begin{array}{ll}\left(n/2\mathrm{,}\left(n-1\right)/2\right)\mathrm{,}\hfill & \text{if}n\text{is}\text{even}\hfill \\ \left(\left(n-1\right)/2\mathrm{,}\left(n-1\right)/2\right)\mathrm{,}\hfill & \text{if}n\text{is}\text{odd}\hfill \end{array}$ (2.2)

The truncated Fourier cosine series can now be written as

$\begin{array}{c}{C}_n\left(x\right)={\displaystyle \sum _{k=0}^n}{f}_k\cos \left(\pi kx\right)\\ ={\displaystyle \sum _{k=0}^{n_1}}{f}_{2k}\cos \left(2\pi kx\right)+{\displaystyle \sum _{k=0}^{n_2}}{f}_{2k+1}\cos \left(\pi \left(2k+1\right)x\right),x\in \left[0,1\right],n\ge 1.\end{array}$ (2.3)

Consider the vector-function

$c\left(\lambda \mathrm{,}x\right)=\left\{-\frac{\pi \csc \left(\frac{\pi \lambda }{2}\right)\cos \left(\pi \lambda \left(x-\frac{1}{2}\right)\right)}{4\lambda }\mathrm{,}-\frac{\pi \sec \left(\frac{\pi \lambda }{2}\right)\sin \left(\pi \lambda \left(x-\frac{1}{2}\right)\right)}{4\lambda }\right\}\mathrm{.}$ (2.4)

It is easy to verify that for a non-integer $\lambda \in ℂ$

$2{\displaystyle {\int }_0^1}c\left(\lambda ,x\right)\cos \left(\pi kx\right)\text{d}x=\left\{\frac{1+{\left(-1\right)}^k}{2\left(k^2-{\lambda }^2\right)},\frac{1+{\left(-1\right)}^k}{2\left(k^2-{\lambda }^2\right)}\right\},k=1,2\dots$ (2.5)

From the Formula (2.1), it follows that the expansion (2.3) is the following representation space $L_2$ as an orthogonal sum of spaces ${\stackrel{1}{L}}_2$ and ${\stackrel{2}{L}}_2$ .

$L_2\left(\mathrm{0,1}\right)={\stackrel{1}{L}}_2\left(\mathrm{0,1}\right)\oplus {\stackrel{2}{L}}_2\left(\mathrm{0,1}\right)\mathrm{.}$ (2.6)

This partition is the starting point of our approach to trigonometric Fourier series.

2.2. Fourier Sine Series

Here

$f\left(x\right)={\displaystyle \sum _{k=1}^{\infty }}{f}_k\sin \left(\pi kx\right)\mathrm{,}x\in \left[\mathrm{0,1}\right]\mathrm{,}$ (2.7)

where

$f\in L_2\left[\mathrm{0,1}\right]\mathrm{,}{f}_s=2{\displaystyle {\int }_0^1}f\left(t\right)\sin \left(\pi kt\right)\text{d}t\mathrm{,}s=\mathrm{1,2,}\cdots$

Now denote for an integer $n\ge 0$

$\left(n_1\mathrm{,}{n}_2\right)=(\begin{array}{ll}\left(n/2-\mathrm{1,}n/2-1\right)\mathrm{,}\hfill & \text{if}n\text{is}\text{even}\hfill \\ \left(\left(n-1\right)/2-\mathrm{1,}\left(n-1\right)/2\right)\mathrm{,}\hfill & \text{if}n\text{is}\text{odd}\hfill \end{array}$ (2.8)

Here the truncated sine Fourier series has the form.

$\begin{array}{c}{S}_n\left(x\right)={\displaystyle \sum _{k=1}^n}{f}_k\sin \left(\pi kx\right)\\ ={\displaystyle \sum _{k=1}^{n_1}}{f}_{2k}\sin \left(2\pi kx\right)+{\displaystyle \sum _{k=0}^{n_2}}{f}_{2k+1}\sin \left(\pi \left(2k+1\right)x\right),x\in \left[0,1\right],n\ge 1.\end{array}$

Consider the vector-function

$s\left(\lambda ,x\right)=\left\{-\frac{1}{4}\pi \csc \left(\frac{\pi \lambda }{2}\right)\sin \left(\pi \lambda \left(x-\frac{1}{2}\right)\right),\frac{1}{4}\pi \sec \left(\frac{\pi \lambda }{2}\right)\cos \left(\pi \lambda \left(x-\frac{1}{2}\right)\right)\right\}$ (2.9)

We see that again

$2{\displaystyle {\int }_0^1}s\left(\lambda ,x\right)\sin \left(\pi kx\right)\text{d}x=\left\{\frac{\left(1+{\left(-1\right)}^k\right)k}{2\left(k^2-{\lambda }^2\right)},\frac{\left(1+{\left(-1\right)}^k\right)k}{2\left(k^2-{\lambda }^2\right)}\right\},k=1,2,\cdots$ (2.10)

Here, we also come to the representation of the form (2.6), in which “sin” appears instead of “cos”.

2.3. Parametric Biorthogonal Systems

Let integers $m\ge 1$ are fixed and is a polynomial.

(2.11)

The coefficients $p_j\in ℂ\mathrm{,}j=\mathrm{0,}\cdots \mathrm{,}m-1$ are now considered arbitrary parameters.

In this section, we will construct biorthogonal systems that form the basis of our approach.

2.3.1. Fourier Cosine Series

Consider the following two unbounded sequences for positive integer s and $r=n_1-m+\mathrm{1,}\cdots \mathrm{,}{n}_1$ in the first sequence and $r=n_2-m+\mathrm{1,}\cdots \mathrm{,}{n}_2$ in the second (see (2.2)).

(2.12)

It’s not hard to see that ${\stackrel{j}{\theta }}_{r\mathrm{,}r}=1$ and for $r\ne s$ ${\stackrel{j}{\theta }}_{r\mathrm{,}s}=0$ , $j=\mathrm{1,2}$ .

Remark 1. For $s=0$ , instead of ${\stackrel{1}{\theta }}_{r\mathrm{,}s}$ , it needs to take ${\stackrel{1}{\theta }}_{r\mathrm{,}s}\mathrm{<mo>/</mo>2}$ . This circumstance is inherent only for this case, and the singularity in the previous formulas for the cosine series at $\lambda =0$ does not disappear (see Remark 4).

For a fixed r and $s\to \infty$ we have ${\stackrel{j}{\theta }}_{r\mathrm{,}s}=O\left(s^{-2}\right)\mathrm{,}j=\mathrm{1,2}$ . Let’s look at the functions

$\begin{array}{l}{\stackrel{1}{\Theta }}_r\left(x\right)=\cos \left(2\pi rx\right)+{\displaystyle \sum _{s=n_1-m}^{\infty }}{\stackrel{1}{\theta }}_{r,s}\cos \left(2\pi sx\right),\\ {\stackrel{2}{\Theta }}_r\left(x\right)=\cos \left(\pi \left(2r+1\right)x\right)+{\displaystyle \sum _{s=n_2-m}^{\infty }}{\stackrel{2}{\theta }}_{r,s}\cos \left(\pi \left(2s+1\right)x\right).\end{array}$ (2.13)

Here, we have two biotogonal and normalized systems, respectively, at ${\stackrel{1}{L}}_2$ and ${\stackrel{2}{L}}_2$ (see (2.6)).

$\left\{\cos \left(2\pi rx\right)\mathrm{,}{\stackrel{1}{\Theta }}_r\left(x\right)\right\}\mathrm{,}\left\{\cos \left(\pi \left(2r+1\right)x\right)\mathrm{,}{\stackrel{2}{\Theta }}_r\left(x\right)\right\}\mathrm{,}x\in \left[\mathrm{0,1}\right]\mathrm{.}$(2.14)

We emphasize that at this stage, this system contains arbitrary parameters $\left\{p_j\right\}$ (see (2.11)).

2.3.2. Fourier Sine Series

Similar to Formula (2.12) (see (2.8))

(2.15)

Similar to the previous case ${\stackrel{j}{\theta }}_{r\mathrm{,}r}=1$ and for $r\ne s$ ${\stackrel{j}{\theta }}_{r\mathrm{,}s}=0$ , $j=\mathrm{1,2}$ . For $s\to \infty$ we have $\stackrel{j}{\theta }=O\left(s^{-1}\right)$ . Here, we will look at the functions.

$\begin{array}{l}{\stackrel{1}{\Theta }}_r\left(x\right)=\sin \left(2\pi rx\right)+{\displaystyle \sum _{s=n_1-m}^{\infty }}{\stackrel{1}{\theta }}_{r,s}\sin \left(2\pi sx\right),\\ {\stackrel{2}{\Theta }}_r\left(x\right)=\sin \left(\pi \left(2r+1\right)x\right)+{\displaystyle \sum _{s=n_2-m}^{\infty }}{\stackrel{2}{\theta }}_{r,s}\sin \left(\pi \left(2s+1\right)x\right)\end{array}$ (2.16)

Here we have two biotogonal and normalized systems for sine, respectively, at ${\stackrel{1}{L}}_2$ and ${\stackrel{2}{L}}_2$ .

$\left\{\sin \left(2\pi rx\right)\mathrm{,}{\stackrel{1}{\Theta }}_r\left(x\right)\right\}\mathrm{,}\left\{\sin \left(\pi \left(2r+1\right)x\right)\mathrm{,}{\stackrel{2}{\Theta }}_r\left(x\right)\right\}\mathrm{,}x\in \left[\mathrm{0,1}\right]\mathrm{.}$

2.3.3. Approximation Formulas for Cosine/Sine Series

When dealing with cosine series, our formula for summation takes the following form (see Section 2.3.1)

$\begin{array}{l}f\left(x\right)\simeq A_{n\mathrm{,}m}\left(x\right)\\ ={\displaystyle \sum _{k=0}^{n_1-m}}{f}_{2k}\cos \left(2\pi kx\right)+{\displaystyle \sum _{k=0}^{n_2-m}}{f}_{2k+1}\cos \left(\pi \left(2k+1\right)x\right)\\ +{\displaystyle \sum _{r=n_1-m+1}^{n_1}}{f}_{2r}{\stackrel{1}{\Theta }}_r\left(x\right)+{\displaystyle \sum _{r=n_2-m+1}^{n_2}}{f}_{2r+1}{\stackrel{2}{\Theta }}_r\left(x\right),x\in \left[0,1\right],n\ge 1.\end{array}$(2.17)

It is easy to see that if we only know the Fourier coefficients $\left\{f_s\right\}\mathrm{,}s\le n$ , the “error of summation” has the following form.

$\begin{array}{l}{\Vert f\left(x\right)-A_{n,m}\left(x\right)\Vert }^2\\ ={\displaystyle \sum _{s=n_1+1}^{\infty }}{\left|f_{2s}-{\displaystyle \sum _{r=n_1-m+1}^{n_1}}{f}_{2r}{\stackrel{1}{\theta }}_{r,s}\right|}^2+{\displaystyle \sum _{s=n_2+1}^{\infty }}{\left|f_{2s+1}-{\displaystyle \sum _{r=n_2-m+1}^{n_2}}{f}_{2r+1}{\stackrel{2}{\theta }}_{r,s}\right|}^2\end{array}$ (2.18)

Similarly, in the case of sine series (see Section 2.3.2)

$\begin{array}{l}f\left(x\right)\simeq B_{n\mathrm{,}m}\left(x\right)\\ ={\displaystyle \sum _{k=0}^{n_1-m}}{f}_{2k}\sin \left(2\pi kx\right)+{\displaystyle \sum _{k=0}^{n_2-m}}{f}_{2k+1}\sin \left(\pi \left(2k+1\right)x\right)\\ +{\displaystyle \sum _{r=n_1-m+1}^{n_1}}{f}_{2r}{\stackrel{1}{\Theta }}_r\left(x\right)+{\displaystyle \sum _{r=n_2-m+1}^{n_2}}{f}_{2r+1}{\stackrel{2}{\Theta }}_r\left(x\right),x\in \left[0,1\right],n\ge 1.\end{array}$(2.19)

$\begin{array}{l}{\Vert f\left(x\right)-B_{n,m}\left(x\right)\Vert }^2\\ ={\displaystyle \sum _{s=n_1+1}^{\infty }}{\left|f_{2s}-{\displaystyle \sum _{r=n_1-m+1}^{n_1}}{f}_{2r}{\stackrel{1}{\theta }}_{r,s}\right|}^2+{\displaystyle \sum _{s=n_2+1}^{\infty }}{\left|f_{2s+1}-{\displaystyle \sum _{r=n_2-m+1}^{n_2}}{f}_{2r+1}{\stackrel{2}{\theta }}_{r,s}\right|}^2\end{array}$(2.20)

3. Summation Algorithms

We proceed from the point of view that the smoother the function that is expanded into one of the series considered here, the more its various coefficients are related to each other. Below, we’d like to present approximation algorithms built based on this spectral approach and the formulas above.

Given the many aspects of the analogy between the sine and cosine series, sometimes we will give all the details only for the cosine series, limiting ourselves to explanations. Sections 3.1 - 3.6 below contain this work’s most important theoretical results, leading to super-fast convergence algorithms for the studied series.

3.1. Basic Linear Equations

3.1.1. The Cosine Series

Let us consider the case of the cosine series. According to the orthogonal expansion (2.6) we have for $f\in {\stackrel{1}{L}}_2\left[\mathrm{0,1}\right]$ , $f=\stackrel{1}{f}+\stackrel{2}{f}$ , $\stackrel{1}{f}\in {\stackrel{1}{L}}_2$ , $\stackrel{2}{f}\in {\stackrel{2}{L}}_2$ .

On the other hand for sine series

$\stackrel{1}{f}\left(x\right)={\displaystyle \sum _{k=0}^{\infty }}{f}_{2k}\cos \left(2\pi kx\right),x\in \left[0,1\right]$

we have (see Section 2.3.3)

${\stackrel{1}{A}}_{n,m}\left(x\right)={\displaystyle \sum _{k=0}^{n_1-m}}{f}_{2k}\cos \left(2\pi kx\right)+{\displaystyle \sum _{r=n_1-m+1}^{n_1}}{f}_{2r}{\stackrel{1}{\Theta }}_r\left(x\right),x\in \left[0,1\right].$

Therefore (see Section 2.3.1)

$\begin{array}{l}\stackrel{1}{f}\left(x\right)-{\stackrel{1}{A}}_{n,m}\left(x\right)\\ ={\displaystyle \sum _{s=0}^{n_1-m}}\left(f_{2s}-{\displaystyle \sum _{r=n_1-m+1}^{n_1}}{f}_{2r}{\stackrel{1}{\theta }}_{r,s}\right)+{\displaystyle \sum _{s=n+1}^{\infty }}\left(f_{2s}-{\displaystyle \sum _{r=n_1-m+1}^{n_1}}{f}_{2r}{\stackrel{1}{\theta }}_{r,s}\right),x\in \left[0,1\right].\end{array}$

From here

$\begin{array}{l}{\Vert \stackrel{1}{f}\left(x\right)-{\stackrel{1}{A}}_{n,m}\left(x\right)\Vert }^2\\ ={\displaystyle \sum _{s=0}^{n_1-m}}{\left|f_{2s}-{\displaystyle \sum _{r=n_1-m+1}^{n_1}}{f}_{2r}{\stackrel{1}{\theta }}_{r,s}\right|}^2+{\displaystyle \sum _{s=n+1}^{\infty }}{\left|f_{2s}-{\displaystyle \sum _{r=n_1-m+1}^{n_1}}{f}_{2r}{\stackrel{1}{\theta }}_{r,s}\right|}^2,x\in \left[0,1\right].\end{array}$(3.1)

The members of the sum at right are unknown to us. We minimize the sum on the left, at least partially (see the beginning of this Section).

Let’s choose a set of m even different coefficients ${\Sigma }_m\subseteq \left\{\mathrm{0,2,}\cdots \mathrm{,2}\left(n_1-m\right)\right\}$ . This is only possible at $m\le n_1/2$ . Our goal is to minimize the sum

${\displaystyle \sum _{s\in {\Sigma }_m}}{\left|f_{2s}-{\displaystyle \sum _{r=n_1-m+1}^{n_1}}{f}_{2r}{\stackrel{1}{\theta }}_{r,s}\right|}^2\to \min .$

in the norm of ${\stackrel{1}{L}}_2\left[\mathrm{0,1}\right]$ concerning coefficients $\left\{p_j\right\}$ of the polynomial (see (2.11)). More details, the following amount is minimized

(3.2)

This leads (see numerators of this amount) to the pseudo-solution of the following linear system

(3.3)

In the ${\stackrel{2}{L}}_2$ we have in corresponding way

${\displaystyle \sum _{s\in {\Sigma }_m}}{\left|f_{2s+1}-{\displaystyle \sum _{r=n_1-m+1}^{n_1}}{f}_{2r+1}{\stackrel{2}{\theta }}_{r,s}\right|}^2\to \min .$

in the norm of ${\stackrel{2}{L}}_2\left[\mathrm{0,1}\right]$ concerning coefficients $\left\{p_j\right\}$ of the polynomial . Here, the corresponding linear system has the form

(3.4)

Note that here ${\Sigma }_m$ consists of m different odd coefficients: ${\Sigma }_m\subseteq \left\{\mathrm{1,3,}\cdots \mathrm{,2}\left(n_1-m\right)+1\right\}$ , $m\le n_2/2$ .

3.1.2. The Sine Series

Similar to the previous case, here we get two linear systems

(3.5)

where ${\Sigma }_m\subseteq \left\{\mathrm{2,4,}\cdots \mathrm{,2}\left(n_1-m\right)\right\}$ , $m\le n_1/2$ .

(3.6)

where ${\Sigma }_m\subseteq \left\{\mathrm{1,3,}\cdots \mathrm{,2}\left(n_2-m\right)+1\right\}$ , $m\le n_2/2$ .

Remark 2. The resulting four equations are essential in our approach (see the beginning of Section 3). Let’s explain some details.

• The determinant of a linear system’s matrix may be zero. To avoid this, one can increase the number of members of the corresponding set ${\Sigma }_m$ to achieve maximum rank for the equation matrix.

• The choice of set ${\Sigma }_m$ affects the algorithm’s results. In our experiments, we used as large numbers as possible (see Section 4.1.2 below).

• We can initially choose different values m₁ and m₂ for ${\stackrel{1}{\theta }}_{r\mathrm{,}s}$ and ${\stackrel{2}{\theta }}_{r\mathrm{,}s}$ correspondingly if desired (see (2.12) in the case of cosine series).

3.2. Definition of Basic Parameters

To minimize the norms mentioned above, linear systems (3.3) - (3.6) must be solved using the method of least squares. When a discrepancy occurs, we record them. Let’s consider the case of a cosine series. We denote the zeros of the

polynomial in (3.3) by $\left\{{\stackrel{1}{\mu }}_j={\left({\stackrel{1}{\lambda }}_{2j}\right)}^2\right\}$ and the zeros of the polynomial in (3.4) by $\left\{{\stackrel{1}{\mu }}_j={\left({\stackrel{2}{\lambda }}_{2j+1}\right)}^2\right\}$ . Note that these polynomials are different. As a result, we have

$\begin{array}{l}{\stackrel{1}{\theta }}_{r,s}=\frac{4r^2-{\left({\stackrel{1}{\lambda }}_{2r}\right)}^2}{4s^2-{\left({\stackrel{1}{\lambda }}_{2r}\right)}^2}{\displaystyle \prod _{\begin{array}{c}p=n_1-m+1\\ p\ne r\end{array}}^{n_1}}\frac{s^2-p^2}{r^2-p^2}\frac{4r^2-{\left({\stackrel{1}{\lambda }}_{2p}\right)}^2}{4s^2-{\left({\stackrel{1}{\lambda }}_{2p}\right)}^2},\\ {\stackrel{2}{\theta }}_{r,s}=\frac{{\left(2r+1\right)}^2-{\left({\stackrel{2}{\lambda }}_{2r+1}\right)}^2}{{\left(2s+1\right)}^2-{\left({\stackrel{2}{\lambda }}_{2r+1}\right)}^2}{\displaystyle \prod _{\begin{array}{c}p=n_2-m+1\\ p\ne r\end{array}}^{n_2}}\frac{{\left(2s+1\right)}^2-{\left(2p+1\right)}^2}{{\left(2r+1\right)}^2-{\left(2p+1\right)}^2}\frac{{\left(2r+1\right)}^2-{\left({\stackrel{2}{\lambda }}_{2p+1}\right)}^2}{{\left(2s+1\right)}^2-{\left({\stackrel{2}{\lambda }}_{2p+1}\right)}^2}.\end{array}$ (3.7)

Remark 3. The notation used here for the zeros of polynomials is convenient because, for example, we see that if ${\stackrel{1}{\lambda }}_{2p}=2p$ at $p\ne r$ , then in ${\stackrel{1}{\theta }}_{r\mathrm{,}s}$ , a cancellation occurs, and instead of m, we have in fact $m-1$ .

3.3. Over-Convergence Generating Spaces

In Section 1.5 we noted that knowledge of a finite number of coefficients of the exponential Fourier series of a function f from the infinite-dimensional space of quasi-polynomials leads to an exact reconstruction of f.

Here, we have four similar spaces. Let’s start with the cosine series (see Section 2.1).

Definition 3. Let $m\ge 1$ be a fixed integer. Consider a system of functions ${\stackrel{1}{U}}_m=\left\{\cos \left(\pi {\lambda }_k\left(x-\frac{1}{2}\right)\right)\right\}$ , ${\lambda }_k\in ℂ$ , $x\in \left[\mathrm{0,1}\right]$ , $k=\mathrm{1,2,}\cdots \mathrm{,}m$ , where $\left\{{\lambda }_k\right\}$ are arbitrary parameters. Consider the linear span ${\stackrel{1}{Q}}_m=span\left\{{\stackrel{1}{U}}_m\right\}$ . We call a function $q\in {\stackrel{1}{Q}}_m$ as a first quasi-cosine of degree at most m.

Similar to Definition 2, it is clear that $q\in {\stackrel{1}{Q}}_m$ if and only if either $q\left(x\right)\equiv 0$ or $q\left(x\right)={\displaystyle {\sum }_k}{P}_{{\beta }_k}\left(x\right)\cos /\sin \left(\pi {\lambda }_{k}x\right)$ , where the polynomials $P_{{\beta }_k}\left(x\right)\cancel{\equiv }0$ are of degree exactly ${\beta }_k$ , and ${\displaystyle {\sum }_k}\left(1+{\beta }_k\right)\le m$ (sin/cos here means: “either sine or cosine”).

Likewise

Definition 4. Let $m\ge 1$ be a fixed integer. Consider a system of functions ${\stackrel{2}{U}}_m=\left\{\sin \left(\pi {\lambda }_k\left(x-\frac{1}{2}\right)\right)\right\}$ , ${\lambda }_k\in ℂ$ , $x\in \left[\mathrm{0,1}\right]$ , $k=\mathrm{1,2,}\cdots \mathrm{,}m$ , where $\left\{{\lambda }_k\right\}$ are arbitrary parameters. Consider the linear span ${\stackrel{2}{U}}_m=span\left\{{\stackrel{2}{V}}_m\right\}$ . We call a function $q\in {\stackrel{2}{Q}}_m$ as a second quasi-cosine of degree at most m.

Similarly, we define the spaces of the first and second-kind quasi-sines. Here the sines and cosines are simply swapped in the same two spaces.

3.4. Explicit Formulas for ${\stackrel{j}{\Theta }}_r$ -Functions

An essential detail of our approach is that the functions ${\stackrel{j}{\Theta }}_r\left(x\right)$ used above have an explicit, elementary form. The evidence for this is similar to the case of the Fourier exponential series (see [25]-[30]).

Let’s focus on the case of space ${\stackrel{1}{Q}}_m$ (see previous section). In Formula (2.12), the right-hand side depends only on s². Let us denote $z=4s^2$ and ${\stackrel{1}{\theta }}_{r\mathrm{,}s}={\stackrel{1}{t}}_{r\mathrm{,}z}$ . The functions ${\stackrel{1}{t}}_{r\mathrm{,}z}$ , considered for $z\in ℂ$ , are rational with ordinary poles at $z={\left({\stackrel{1}{\lambda }}_{2p}\right)}^2$ . Let $U\subset ℂ$ be a bounded simply connected open subset containing all points $\left\{{\stackrel{1}{\lambda }}_{2p}\right\}$ , with the positively oriented simple boundary curve $\gamma =\partial U$ .

We have ${\stackrel{1}{t}}_{r\mathrm{,}z}=O\left(1/z\right)$ , $z\to \infty$ , therefore, according to Cauchy’s residue theorem

$0=\frac{1}{2\pi i}{\displaystyle {\int }_{\gamma }}\frac{{\stackrel{1}{t}}_{r\mathrm{,}t}}{t-s^2}\text{d}t={\stackrel{1}{t}}_{r\mathrm{,}{s}^2}+{\displaystyle \sum _{\forall p}}{\text{Res}}_{z={\left({\stackrel{1}{\lambda }}_{2p}\right)}^2}\left(\frac{{\stackrel{1}{t}}_{r\mathrm{,}z}}{z-s^2}\right)\mathrm{,}s\in U\backslash \left\{{\stackrel{1}{\lambda }}_{2p}\right\}\mathrm{.}$

So

${\stackrel{1}{t}}_{r\mathrm{,}s}=-{\displaystyle \sum _{\forall p}}{\text{Res}}_{z={\left({\stackrel{1}{\lambda }}_{2p}\right)}^2}\left(\frac{{\stackrel{1}{t}}_{r\mathrm{,}z}}{z-s^2}\right)\mathrm{,}s\in U\backslash \left\{{\stackrel{1}{\lambda }}_{2p}\right\}\mathrm{.}$(3.8)

These residues can be explicitly calculated. To avoid clattering by detailed formulas in the presentation, we present corresponding results only for cosine series in ${\stackrel{1}{L}}_2$ , when all parameters $\left\{{\stackrel{1}{\lambda }}_{2p}\right\}$ are different.

Theorem 2. When for $p\ne q$ , ${\stackrel{1}{\lambda }}_{2p}\ne {\stackrel{1}{\lambda }}_{2q}$ , the system $\left\{{\stackrel{1}{\Theta }}_r\right\}$ is expressed directly in the form

${\stackrel{1}{\Theta }}_r\left(x\right)={\displaystyle \sum _{k=n_1-m+1}^{n_1}}{\stackrel{1}{c}}_{r\mathrm{,}k}\cos \left(\pi {\stackrel{1}{\lambda }}_{2k}\left(x-1/2\right)\right)\mathrm{,}x\in \left[\mathrm{0,1}\right]\mathrm{,}$(3.9)

where $r=n_1-m+\mathrm{1,}\cdots \mathrm{,}{n}_1$ and³

$\begin{array}{l}{\stackrel{1}{c}}_{r,k}=\frac{\pi \left(r+{\stackrel{1}{\lambda }}_{2k}/2\right)}{4{\stackrel{1}{\lambda }}_{2k}\text{sinc}\left(\pi \left(r-{\stackrel{1}{\lambda }}_{2k}/2\right)\right)}{\displaystyle \prod _{\begin{array}{c}p=n_1-m+1\\ p\ne k\end{array}}^{n_1}}\frac{4r^2-{\left({\stackrel{1}{\lambda }}_{2k}\right)}^2}{{\left({\stackrel{1}{\lambda }}_{2k}\right)}^2-{\left({\stackrel{1}{\lambda }}_{2p}\right)}^2}\\ {\displaystyle \prod _{\begin{array}{c}q=n_1-m+1\\ q\ne r\end{array}}^{n_1}}\frac{{\left({\stackrel{1}{\lambda }}_{2k}\right)}^2-4q^2}{4r^2-4q^2},r,k=n_1-m+1,\cdots ,n_1.\end{array}$

The proof is similar to the case of the exponential Fourier series (see [25] [27]). Similar formulas are valid in the other three cases (for ${\stackrel{2}{\Theta }}_r\left(x\right)$ in the cosine series and ${\stackrel{j}{\Theta }}_r\left(x\right)\mathrm{,}\left(j=\mathrm{1,2}\right)$ in the case of the sine series.

Remark 4. In the above three cases, all singularities are eliminated using the sinc-function. However, Formula (3.9) shows that at ${\stackrel{1}{\lambda }}_{2k}=0$ in the cosine series, the singularity is preserved. For a solution to this problem numerically, please look at Section 4.1.3 below.

3.5. Adaptive Summation Formulas

Let’s denote (see above Section)

${\stackrel{1}{a}}_k={\displaystyle \sum _{r=-m+n_1+1}^{n_1}}{f}_{2r}{\stackrel{1}{c}}_{r,k},{\stackrel{2}{a}}_k={\displaystyle \sum _{r=-m+n_2+1}^{n_2}}{f}_{2r+1}{\stackrel{2}{c}}_{r,k}.$(3.10)

Consider functions (see Section 2.1 and Theorem 2)

$\begin{array}{l}{\stackrel{1}{O}}_{n,m}\left(x\right)={\displaystyle \sum _{k=-m+n_1+1}^{n_1}}{\stackrel{1}{a}}_k\cos \left(\pi {\stackrel{1}{\lambda }}_{2k}\left(x-\frac{1}{2}\right)\right),x\in \left[0,1\right]\\ {\stackrel{2}{O}}_{n,m}\left(x\right)={\displaystyle \sum _{k=-m+n_1+1}^{n_1}}{\stackrel{2}{a}}_k\cos \left(\pi {\stackrel{2}{\lambda }}_{2k+1}\left(x-\frac{1}{2}\right)\right),x\in \left[0,1\right].\end{array}$(3.11)

Let’s denote now

$\begin{array}{l}{\stackrel{1}{A}}_{n\mathrm{,}m}\left(x\right)={\displaystyle \sum _{s=0}^{n_1-m}}\left(f_{2s}-{\displaystyle \sum _{r=-m+n_1+1}^{n_1}}{f}_{2r}{\stackrel{1}{\theta }}_{r\mathrm{,}s}\right)\cos \left(2\pi sx\right)+{\stackrel{1}{O}}_{n\mathrm{,}m}\left(x\right)\mathrm{,}x\in \left[\mathrm{0,1}\right]\\ {\stackrel{2}{A}}_{n\mathrm{,}m}\left(x\right)={\displaystyle \sum _{s=0}^{n_2-m}}\left(f_{2s+1}-{\displaystyle \sum _{r=-m+n_2+1}^{n_2}}{f}_{2r+1}{\stackrel{2}{\theta }}_{r\mathrm{,}s}\right)\cos \left(\left(2s+1\right)\pi x\right)+{\stackrel{2}{O}}_{n\mathrm{,}m}\left(x\right)\mathrm{,}x\in \left[\mathrm{0,1}\right]\mathrm{.}\end{array}$(3.12)

The function $f\left(x\right)\in L_2\left[\mathrm{0,1}\right]$ is approximated by the formula

$f\left(x\right)\simeq A_{n\mathrm{,}m}\left(x\right)={\stackrel{1}{A}}_{n\mathrm{,}m}\left(x\right)+{\stackrel{2}{A}}_{n\mathrm{,}m}\left(x\right)\mathrm{,}x\in \left[\mathrm{0,1}\right]$(3.13)

Similar formulas were applied to the case of the sine series. This approximations was, in fact, constructed according to Krylov’s scheme (see Section 1.2). We call the algorithms corresponding to the above summation formulas adaptive since the terms ${\stackrel{j}{O}}_{n\mathrm{,}m}\left(x\right)\mathrm{,}j=\mathrm{1,2}$ appearing as the function g from (1.3) depends on the function $f\left(x\right)$ (“adapt” to it). Unlike the KE-method, which can only use a few g-functions (Section 1.3), here a broadest class of trigonometric functions that depend on parameters from $ℂ$ are used (see Section 3.3).

Remark 5. When solving linear equations in Section 3.1, sometimes the determinant of the matrix turns out to be zero. (see Remark 2 and beginning of Section 3.2). Formulas (3.12) take into account the resulting discrepancies. Above in (3.12), this is the next part of the sum

$d_1\left(x\right)={\displaystyle \sum _{s\in {\Sigma }_m}}\left(f_{2s}-{\displaystyle \sum _{r=-m+n_1+1}^{n_1}}{f}_{2r}{\stackrel{1}{\theta }}_{r,s}\right)\cos \left(2\pi sx\right),x\in \left[0,1\right]$

which is equal to zero in $\left[\mathrm{0,1}\right]$ if the determinant differs from zero.

3.6. The Over-Convergence Phenomenon

This main theoretical result of this work is valid for all four proposed algorithms. Let’s look at the case of cosine series with even coefficients (see ${\stackrel{1}{A}}_{n\mathrm{,}m}$ in (3.12)). Denote $D_m=\left\{2\left(n_1-m+1\right)\mathrm{,}\cdots \mathrm{,2}{n}_1\right\}$ and take ${\Sigma }_m$ from (3.5).

Theorem 3 (The phenomenon of the over-convergence). Let $f\in {\stackrel{1}{Q}}_m$ (see Definition 3), and the Fourier coefficients $\left\{f_{2s}\right\}$ , $s\in D_m\cup {\Sigma }_m$ , of the function f be given. Denote by Λ the set of all integer parameters in $\left\{{\stackrel{1}{\lambda }}_{2j}\right\}$ (see Section

3.2). In order for the ${\stackrel{1}{A}}_{n\mathrm{,}m}\left(x\right)$ to be exact (that is, $f\left(x\right)\equiv {\stackrel{1}{A}}_{n\mathrm{,}m}\left(x\right)$ , $x\in \left[\mathrm{0,1}\right]$ ), it is necessary and sufficient that $\Lambda \subseteq D_m\cup {\Sigma }_m$ .

If we keep Formula (3.8) in mind, then the proof is quite similar to that given in [25].

Definition 5. The algorithms underlying the summation formulas ${\stackrel{1}{A}}_{n\mathrm{,}m}\left(x\right)$ and ${\stackrel{2}{A}}_{n\mathrm{,}m}\left(x\right)$ (see (3.12)) will be denoted by ${ℭ}_1$ and ${ℭ}_2$ , respectively. We denote the corresponding two algorithms for Fourier sine series by ${\mathfrak{S}}_1$ and ${\mathfrak{S}}_2$ .

Accordingly, we denote the parallel algorithm for the cosine Fourier series by $ℭ$ and the sine series by $\mathfrak{S}$ .

3.7. Fourier Ordinary Series Revisited

3.7.1. Adaptive Method for Fourier Exponential Series

First, we’d like to present briefly a scheme for accelerating the convergence of work [29]. Let $m\ge 1$ be an integer and $D_m$ be a set of m different integers and , $z\in ℂ$ , be a polynomial

(3.14)

where $\left\{{\lambda }_k\right\}$ , $k\in D_m$ , be the set of its zeros. Consider the sequence

(3.15)

(see Section 2.3.2).

$T_r\left(x\right)=\exp \left(i\pi rx\right)+{\displaystyle \sum _{s\notin D_m}}{t}_{r\mathrm{,}s}\exp \left(i\pi sx\right)\mathrm{,}r\in D_m$ (3.16)

The space of quasi-polynomials $Q_m$ is defined (see Definition 2 above), and an explicit form of the functions $T_r\in Q_m$ is given.

The basic system of linear equations regarding the definition of $\left\{a_k\right\}$ (see (3.14)) has the following form

(3.17)

where ${\tilde{D}}_m$ is a set of other m different integers (${\tilde{D}}_m\cap D_m=\varnothing$ ).

The next step is finding the zeros of polynomial . Finally, function $f\left(x\right)$ is approximated by the summation formula

$f\simeq F_{n\mathrm{,}m}\left(x\right)={\displaystyle \sum _{r\in D_m}}{f}_r{T}_r\left(x\right)+{\displaystyle \sum _{s\in D_n\backslash D_m}}\left(f_s-{\displaystyle \sum _{r\in D_m}}{f}_r{t}_{r\mathrm{,}s}\right)\exp \left(i\pi sx\right)\mathrm{.}$ (3.18)

with the error

$R_{n\mathrm{,}m}\left(x\right)\stackrel{def}{=}f\left(x\right)-F_{n\mathrm{,}m}\left(x\right)\mathrm{.}$ (3.19)

This algorithm is called adaptive since equation (3.17) “adapts” solutions $\left\{a_k\right\}$ to the behavior of coefficients $\left\{f_s\right\}$ . It can be noted that this scheme differs from our previous cases only in notations.

3.7.2. About Historically Initial Form of Fourier Series

Joseph Fourier used the trigonometric form of the Fourier series, although Euler’s formula for the exponential was known long before this. Solving a practical problem, instead of approximation (1.1), he studied the formula

$\begin{array}{l}f\left(x\right)\simeq {\displaystyle \sum _{k=-n}^n}\left(a_k\cos \pi kx+b_k\sin \pi kx\right)\mathrm{,}x\in \left[-\mathrm{1,1}\right]\mathrm{.}\\ a_0=\frac{1}{2}{\displaystyle {\int }_{-1}^1}f\left(t\right)\text{d}t,a_k={\displaystyle {\int }_{-1}^1}f\left(t\right)\cos \pi kt\text{d}t,b_k={\displaystyle {\int }_{-1}^1}f\left(t\right)\sin \pi kt\text{d}t,\end{array}$ (3.20)

where f is a real-valued function.

Do the acceleration algorithms proposed above work now similarly to the algorithm for (1.1)? The exact answer is not so trivial. First, we show that the partial sum (1.1) is not generally equivalent to the partial sum of the form (3.20). For example, if in (1.1) we go to trigonometric functions, then for even n, the summation over sines is carried out in the range from 1 to $\left(n-1\right)$ , while over cosines - from 0 to n. In addition, even coefficients for sines and odd - for cosines are equal to zero. Therefore, we can apply algorithm ${ℭ}_1$ to the cosine sum and to the sine sum - ${\mathfrak{S}}_1$ (see Definition 5). Two algorithms are used here in parallel (see Sections 3.5.1 and 3.5.2).

As for the sum (3.20), we can consider the sums of cosines and sines separately on $\left[\mathrm{0,1}\right]$ and continue each of them to $\left[-\mathrm{1,0}\right)$ , even and odd. All four of our algorithms work here in parallel. Thus, the numerical complexity of this algorithm is twice more at the same n.

Recall that the adaptive algorithm for the Fourier exponential series was designated $\mathfrak{A}$ (see [25] and Section 1.5 above). We’d like to introduce notations for the algorithms mentioned above.

Definition 6. For the algorithm of article [25], we will leave the designation $\mathfrak{A}$ . We denote the parallel algorithm for the trigonometric form of the exponential Fourier series by $ℭ{\mathfrak{S}}_{12}$ and the parallel algorithm for (3.20) - by $ℭ\mathfrak{S}$ .

4. The Numerical Experiment

4.1. Implementation Remarks

4.1.1. Testing Functions

Numerical experiments were carried out with the following functions.

$f_1\left(x\right)=(\begin{array}{ll}\mathrm{0,}\hfill & 0\le x\mathrm{<mo><</mo>}1/3\mathrm{<mo>;</mo>}\hfill \\ \frac{29}{20}{\sin }^5\left(3x-1\right)\mathrm{,}\hfill & 1/3\le x\le \mathrm{1<mo>;</mo>}\hfill \end{array}$ (4.1)

$f_2\left(x\right)=\frac{63}{89}{\text{e}}^{2i{\left(x+\frac{2}{3}\right)}^2}\mathrm{,}{f}_3\left(x\right)=\frac{\cos \left(x\right)}{\left(2/3+i\right)-x}\mathrm{,}x\in \left[\mathrm{0,1}\right]$ (4.2)

Here $f_1\left(x\right)\in C^5$ and $f_1{\left(x\right)}^{\left(5\right)}$ is a piece-wise smooth function with two jumps: at $x=1/3$ and $x=1$ (see (1.1)). The other two functions are analytic in a neighborhood of segment [0, 1]. Function $f_2\left(x\right)$ is entire, but not exponential type. Function $f_3\left(z\right)$ is rational and has only one pole at $z=3/2$ . Functions $f_1\left(x\right)$ and $f_2\left(x\right)$ were tested in [29] (on $x\in \left[-\mathrm{1,1}\right]$ with the linear change of variable).

4.1.2. Used Algorithms

There are usually different ways to choose each of the four ${\Sigma }_m$ -sets of coefficients to solve the corresponding equations (see Section 3.1). We used as large numbers as possible. For example, for algorithm ${ℭ}_1$ we have ${\Sigma }_m=\left\{2\left(n_1-2m+1\right)\mathrm{,}\cdots \mathrm{,2}\left(n_1-m\right)\right\}$ .

The parallel algorithms defined above have fewer possibilities for increasing m with the same amount of coefficients used as algorithm $\mathfrak{A}$ . The main algorithms we are studying, $ℭ$ and $\mathfrak{S}$ , have better opportunities in this regard $n_1=n_2$ , since in the best case $m=\min \left(n_1\mathrm{,}{n}_2\right)$ .

Definition 7. Let $m\ge 1$ be fixed. We call algorithms $ℭ$ and $\mathfrak{S}$ optimal, if n is in (2.2) and (2.8), respectively, odd and even, and the minimum number of Fourier coefficients is used.

4.1.3. Calculation Details

Our codes are written in Wolfram Mathematica language (see [31]). The calculations were performed on a standard personal computer. Its were carried out using from 16 to 128 decimal places (from 16 to 128 working precision by Wolfram Mathematica terminology), depending on the value of n.

Using the symbolic capabilities of the system Wolfram Mathematica, we chose the tested functions such that the corresponding Fourier coefficients can be representable in an exact form using special functions. When choosing an algorithm’s working precision, we translate in advance the exact (symbolic) values of the Fourier coefficients into numerical values of the same accuracy.

It is very helpful to use additional processing for the data produced by the algorithm. For example, when applying a working precision equal to w, we consider that two values of parameters ${\stackrel{1}{\lambda }}_{2k}$ and ${\stackrel{1}{\lambda }}_{2p}$ are equal to each other if $\left|{\stackrel{1}{\lambda }}_{2k}-{\stackrel{1}{\lambda }}_{2p}\right|\le {10}^{-w/3}$ . Such a seemingly rigid approach has justified itself in our

experiments, preventing calculations from being interrupted due to the appearance of the uncertainty of type 0/0. The same must be kept in mind in the case where the function f is real and the complex conjugate zeros of the polynomial are calculated with shifts (see Section 2.3), which interferes with the reality of the approximation.

As for the problem with singularity in the cosine series (see Remarks 1 and 4), one can “correct” the value ${\stackrel{1}{\lambda }}_{2k}=0$ to ${\stackrel{1}{\lambda }}_{2k}={10}^{-w/3}$ here, too. The values of w

we use are so large (compared to the expected approximation error) that the algorithm will not notice such “manipulations”.

With integer parameters, the implementation of shortcuts described in Remark 2 can be difficult numerically. Our simple method is quite effective here: designate such parameters symbolically, and only at the algorithm’s output do we use their limits.

4.1.4. Tables Information

Table 1 uses the following notation for partial sums of classical Fourier series: FTrig—Formula (3.20); FExp—(1.1); FCosine—(2.1); FSine—(2.7);

In the Tables below, the following designations are used: Alg.—algorithms; $\mathcal{N}$ —number of Fourier coefficients used.

The main algorithms’ names are used from Definitions 5, 6 and 7.

Relative L₂-error for an approximation $\tilde{f}$ of function f is calculated by the formula

$\Vert f-\tilde{f}\Vert =(\begin{array}{ll}\frac{{\Vert f-\tilde{f}\Vert }_2}{{\Vert f\Vert }_2}\mathrm{,}\hfill & {\Vert f\Vert }_2\ne \mathrm{0<mo>;</mo>}\hfill \\ \mathrm{0,}\hfill & {\Vert f\Vert }_2=0.\hfill \end{array}$ (4.3)

4.2. Numerical Results

Table 1. L₂ errors of function $f_3\left(x\right)$ using classical partial sums of Fourier series.

Alg. $\downarrow$	n = 20	n = 40	n = 60	n = 80	n = 100	n = 120
FTrig	9.5e−2	4.67e−2	4.1e−2	3.7e−2	3.4e−2	3.1e−2
FExp	9.5e−2	6.75e−2	5.52e−2	4.79e−2	4.28e−2	3.91e−2
FSine	1.83e−1	1.29e−1	1.06e−1	9.18e−2	8.21e−2	7.5e−2
FCosine	2.13e−2	7.62e−3	4.17e−3	2.72e−3	1.95e−3	1.48e−3

Table 2. L₂ errors of function $f_2\left(x\right)$ using optimal adaptive algorithms at practically equal $\mathcal{N}$ with $\mathfrak{A}$ .

Alg.	m = 4	m = 8	m = 12	m = 16	m = 20	m = 24	m = 28	m = 32
$\mathfrak{A}$	1e−3	1.2e−8	2.4e−14	1.4e−20	3e−27	3e−34	1.6e−41	4.7e−49
$ℭ$	5.8e−4	1.1e−8	3.2e−14	2.7e−20	8.1e−27	1.1e−33	8.2e−41	3.4e−48
$\mathfrak{S}$	1.2e−2	6.2e−7	3.4e−12	4.3e−18	1.9e−24	3.5e−31	3.2e−38	1.5e−36

Table 3. L₂ errors of function $f_3\left(x\right)$ using optimal adaptive algorithms at practically equal $\mathcal{N}$ with $\mathfrak{A}$ .

Alg.	m = 4	m = 8	m = 12	m = 16	m = 20	m = 24	m = 28	m = 32
$\mathfrak{A}$	2.1e−5	9.2e−11	2.7e−17	2.9e−23	4.5e−29	1.3e−35	3.4e−41	1e−47
$ℭ$	2.1e−5	2e−10	3.1e−15	6e−20	1.3e−24	2.9e−29	6.6e−34	1.5e−38
$\mathfrak{S}$	1.6e−2	1.1e−8	1.8e−13	3.6e−18	9.1e−23	2.5e−27	6.2e−32	1.5e−36

Table 4. L₂-errors for function $f_1\left(x\right)$ at $\mathcal{N}=const.$ and small m.

Alg. $\downarrow$	$\mathcal{N}\downarrow$	m = 1	m = 2	m = 3	m = 4	m = 5	m = 6
$\mathfrak{A}$	25	3.9e−4	2.8e−6	4.2e−5	2.2e−6	3.3e−5	6.8e−6
$ℭ$	26	3.2e−6	6.6e−6	9.9e−7	1.6e−6	5.7e−7	9.6e−7
$\mathfrak{S}$	26	1.9e−5	1.4e−5	1.2e−5	7e−6	7.5e−6	7.1e−6

Table 5. L₂-errors for function $f_2\left(x\right)$ at $\mathcal{N}=const.$ and small m.

Alg. $\downarrow$	$\mathcal{N}\downarrow$	m = 1	m = 2	m = 3	m = 4	m = 5	m = 6
$\mathfrak{A}$	25	2.4e−3	6.6e−6	5.1e−7	5.2e−9	6.8e−10	5.2e−11
$ℭ$	26	1.5e−6	1.5e−9	1.9e−11	6e−13	7.4e−14	3.2e−14
$\mathfrak{S}$	26	7.4e−5	8.5e−8	8.3e−10	3.8e−11	6.7e−12	3.4e−12

4.3. A Couple of Notes after Experiments

Table 1 shows that only the partial sum of the cosine series shows a relatively noticeable practical approximation. In the remaining tables, we demonstrated the operation of our three main adaptive algorithms.

Table 2 and Table 3 show the operation of the algorithms over a relatively “long distance”. We see a clear advantage of $\mathfrak{A}$ in approximation speed, but we can still increase the parameters up to $m=16$ . At a “short distance,” algorithm $ℭ$ retains the lead, and $\mathfrak{S}$ overtakes $\mathfrak{A}$ (see Table 5). As for the approximation of the function of limited smoothness $f_1\left(x\right)$ , the maximum accuracy of the approximation is achieved here already at $m=3$ (Table 4).

We did not include algorithms in the tables that contain combinations of cosine and sine series elements (see Definition 6). Their results are naturally similar to those of the sine series. However, we think it’s necessary to note that such adaptive parallel algorithms can be greatly improved by changing the traditional concept of a partial sum. For example, the summation in the Formula (3.20) can be considered over cosines with smaller n and over sines—with bigger ones without changing the total number of coefficients (here equal to $1+2n$ ).

5. Conclusions

5.1. Summary

The phenomenon of super-convergence for the Fourier series, discovered in work [25], made it possible to obtain fast, high-accuracy spectral algorithms using a parametric biorthogonal system. The latest, modified version of used method was presented and numerically tested in the works [29] and [30].

This work applies this approach to the cosine and sine Fourier series. A new circumstance arose when it turned out that the direct application of the method did not work, and it was necessary first, divide the spectrum into two parts. The algorithms turned out to be parallel. From a computational point of view, this resulted in the fact that instead of solving a linear system plus finding the roots of a polynomial, the order of the two equations and the powers of the two polynomials were halved. Already, the first of these steps reduces the complexity of the algorithm by fourfold. The stage of finding the zeros of the polynomial will also be significantly more accessible (see Section 3.2). For example, if in the case of the Fourier series (algorithm $\mathfrak{A}$ ), the zeros are determined exactly at $m\le 3$ , then here—at $m\le 6$ . However, the “convergence acceleration” by algorithm $\mathfrak{A}$ is much more impressive.

The algorithms based on the phenomenon of over-convergence have at least one drawback: they require excellent smoothness for the approximated function. Apparently, under this condition, the Gibbs phenomenon is overcome here in the best possible way.

Let us present the main results of the work.

• The convergence acceleration method based on a parametric biorthogonal system is theoretically justified and numerically implemented for the trigonometric Fourier series. It is precisely these series and their discrete analogues that are most often used in modern applications: from signal/image processing and computer tomography to physics and astronomy.

• An important aspect here is the detection of the phenomenon of over-convergence (see Section 3.6) which can equally be considered the basis of our method. The potential of this amazing phenomenon is far from exhausted.

• Impressive numerical experiment results (see Section 4), obtained using standard personal computer systems, demonstrate the reliability of the proposed algorithms.

• There is every reason to believe that the developed approach can be successfully applied to many more general problems (see below).

5.2. Future Challenges

Some predictions about future generalizations of algorithm $\mathfrak{A}$ were presented in the papers [29] and [30]. The part of this can be repeated here without change. This applies, in particular, to multidimensional Fourier series and to the case when the approximated function is piecewise smooth.

Let’s look only at some new aspects.

5.2.1. Spectrum Splitting Algorithms for Sturm-Liouville Problem

Direct and inverse Sturm-Liouville problems are seen as the closest “consumers” of our adaptive algorithms. Firstly, it would be interesting to discover the phenomenon of super-convergence here. Secondly (if this is confirmed), it would be natural to apply the corresponding adaptive algorithms to the numerical solution of inverse problems with a smooth potential.

Over the past decades, the works of Vladislav V. Kravchenko and a number of his co-authors have developed effective numerical algorithms for solutions to such problems (see, for example, [32]-[34]). Our approaches would be a natural complement to these studies.

5.2.2. On Boundary Value Problems for General ODE

The most significant area for future research is the following self-adjoint boundary value problem

$Ly=\lambda y,Ly={\displaystyle \sum _{k=0}^n}{\left(p_k\left(x\right)y^{\left(k\right)}\left(x\right)\right)}^{\left(k\right)},\lambda \in ℂ,x\in \left[0,1\right]$ (5.1)

where $n\ge 1$ , $q\gg 1$ , $p_n\left(x\right)=1$ , $p_k\left(x\right)\in C^{k+nq}\left[\mathrm{0,1}\right]$ , with regular boundary conditions $\left\{U_{\mu }=0\right\}$ , $\mu =\mathrm{1,2,}\cdots \mathrm{,2}n$ (for details, see monograph [35]).

Unfortunately, the conclusions in [29] and [30] were too optimistic about the possibilities of the corresponding generalizations.

Possibility to use even greater parallelization at $n\ge 2$ and the asymptotic estimates for the problem (5.1) (see [35]) allows us to see the possibility of detecting the phenomenon of over-convergence here, too.

Acknowledgements

I want to thank an anonymous referee whose objections and comments have significantly improved the exposition of this work.

NOTES

¹Historically more accurate: “Wilbraham-Gibbs phenomenon”.

²The jump at point $a_m=1$ is considered to be $A_{m\mathrm{,}k}\left(f\right)=f^{\left(k\right)}\left(1-0\right)-f^{\left(k\right)}\left(-1+0\right)$ .

³$\text{sinc}\left(t\right)=\sin \left(t\right)/t$ , $\text{sinc}\left(0\right)=1$ .