APPROXIMATION OF TIME FUNCTIONS BY THE INTERPOLATION METHOD

V. M. OSIPOV

(Presented by the scientific-methodological seminar of the Department of Engineering and Computational Mathematics).

The problem of approximating time functions, or other functions defined on a semi-infinite interval, often arises in various branches of modern engineering.

If a function is given analytically, then the approximation problem is best solved by expanding it into a finite series with respect to one or another system of orthogonal functions. In this case, as is known, the approximation error will be minimal in the metric of the corresponding Hilbert space. Among all systems of this kind, the exponential polynomials (“e”-polynomials) of Chebyshev of the first kind \(T_k^*(t)\) \((k=0, 1, 2, \ldots)\) and the “e”-functions of Chebyshev of the third kind \(S_k(t)\) \((k=1, 2, \ldots)\) possess the best approximating capabilities [5]. However, if the Fourier coefficients of the expansion in the functions \(S_k(t)\) \((k=1, 2, \ldots)\) can be obtained by means of exponential moments for a very broad class of time functions [6], then the expansion in the “e”-polynomials \(T_k^*(t)\) \((k=0, 1, 2, \ldots)\), as a rule, cannot be found [5]. Meanwhile, approximate values of the coefficients of both expansions, differing only slightly from the Fourier coefficients, can be obtained as a result of an interpolation process in which the interpolation nodes are the zeros of the “e”-polynomial \(T_n^*(t)\) or of the function \(S_{n+1}(t)\). Such a way of determining the expansion coefficients is the simplest and entirely natural when the time function to be expanded is given in the form of some empirical dependence, i.e., graphically.

A remarkable property of ordinary orthogonal polynomials defined on a finite interval—namely, that they are orthogonal not only in the sense of integration but also in the sense of summation over the interpolation nodes [4]—extends without any changes to the Chebyshev “e”-functions \(T_k^*(t)\) and \(S_k(t)\), which substantially simplifies the computational procedure.

In this connection, interpolation problems with the aid of orthogonal functions are conveniently considered from the standpoint of the simplest concepts of functional analysis.

Various interpolation schemes using the “e”-functions \(T_k^*(t)\) and \(S_k(t)\) \((k=0, 1, 2, \ldots, n)\) are possible; to each of these schemes there corresponds a unique orthogonal basis generating some subspace of the “\(n\)”-dimensional Euclidean space \(R_n\). Some of these schemes are considered below.

1. Orthogonal bases of certain subspaces

a) The subspace \(R_n^{\mathrm{tt}}\)

Let \(T_0^*(t), T_1^*(t), \ldots, T_{n-1}^*(t)\) be the first “\(n\)” Chebyshev “\(e\)”-polynomials of the first kind [5]. Suppose that \(t_i^{\mathrm{t}}\) \((i=1,2,\ldots,n)\) are the zeros of the “\(e\)”-polynomial \(T_n^*(t)\). We form a system of \(n\) vectors \(T_k\) \((k=0,1,\ldots,n-1)\). The \(n\) coordinates of each vector of this system are defined as the values of the “\(e\)”-polynomial \(T_k^*(t)\) at the zeros of the “\(e\)”-polynomial \(T_n^*(t)\), i.e.,
\[ T_k=\{T_k^*(t_1^{\mathrm{t}}),\ T_k^*(t_2^{\mathrm{t}}),\ \ldots,\ T_k^*(t_n^{\mathrm{t}})\}. \tag{1-1} \]

The constructed system of vectors forms an orthogonal basis of a certain subspace \(R_n^{\mathrm{tt}}\) of the Euclidean space \(R_n\) [1], since [Appendix 1]
\[ (T_k,T_m)=0 \quad k\ne m \]
\[ \|T_k\|^2=(T_k,T_k)=\frac{n}{2}. \tag{1-2} \]

Any vector \(f^{\mathrm{t}}\in R_n^{\mathrm{tt}}\) can be represented in the form [1]
\[ f^{\mathrm{t}}=\sum_{k=0}^{n-1} b_k^{\mathrm{t}} T_k. \tag{1-3} \]

In particular, suppose that a certain bounded function \(f(t)\) \((0\le t\le \infty)\) is given and satisfies the condition [5]
\[ \int_0^\infty \frac{e^{-\frac{at}{2}}[f(t)]^2}{\sqrt{1-e^{-a}}}\,dt<\infty, \]
then it can be expanded in a convergent series in the “\(e\)”-polynomials \(T_k^*(t)\)
\[ f(t)=\frac{b_0}{2}+\sum_{k=1}^{\infty} b_k T_k^*(t). \]

Consider a segment of the series
\[ f(t)\simeq \frac{b_0}{2}+\sum_{k=1}^{n-1} b_k T_k^*(t). \tag{1-4} \]

Set successively \(t=t_i^{\mathrm{t}}\) \((i=1,2,\ldots,n)\); then we obtain a system of \(n\) linear equations, which in vector form can be written as follows:
\[ f^{\mathrm{t}}=\sum_{k=0}^{n-1} b_k^{\mathrm{t}} T_k, \tag{1-5} \]
where the vector \(f^{\mathrm{t}}\) has the form
\[ f^{\mathrm{t}}=\{f(t_1^{\mathrm{t}}),\ f(t_2^{\mathrm{t}}),\ \ldots,\ f(t_n^{\mathrm{t}})\}, \tag{1-6} \]
\(b_k^{\mathrm{t}}\) \((k=0,1,\ldots,n-1)\) are scalar quantities, differing somewhat from the coefficients \(b_k\) \((k=0,1,\ldots,n-1)\) in the expansion (1-4); \(T_k\) \((k=0,1,\ldots,n-1)\) is the system of vectors defined by formula (1-1) and forming an orthogonal basis of the subspace \(R_n^{\mathrm{tt}}\), with the vector \(T_0\) being taken with weight \(1/2\).

Equality (1-5) coincides with (1-3). It means that there exists an operator \(\tilde I_n^{\mathrm{tt}}\) effecting an isometric mapping of the subspace \(R_n^{\mathrm{tt}}\) onto the \(n\)-dimensional Euclidean space \(R_n\), i.e.,
\[ \tilde I_n^{\mathrm{tt}}\cdot f^{\mathrm{t}}=b^{\mathrm{t}}, \tag{1-7} \]
where
\[ b^{\mathrm{t}}=(b_0^{\mathrm{t}},\ b_1^{\mathrm{t}},\ \ldots,\ b_{n-1}^{\mathrm{t}}). \tag{1-8} \]

Let us clarify the structure of the operator \(\widetilde{I}_n^{\tau\tau}\). Multiply (1-5) scalarly by the vector \(T_m\); then, by virtue of (1-2), we shall have

\[ b_m^\tau=\frac{(f^\tau,T_m)}{\|T_m\|^2} =\frac{2}{n}\sum_{i=1}^{n} f(t_i^\tau)T_m^*(t_i^\tau)(m=0,1,2\ldots n-1). \tag{1-9} \]

The set of these equations can be written in the form (1-7). Consequently, the operator \(\widetilde{I}_n^{\tau\tau}\) is determined by an \(n\times n\) matrix with elements

\[ a_{ki}=\frac{T_k^*(t_i^\tau)}{\|T_k\|^2} =\frac{2}{n}T_k^*(t_i^\tau)(k=0,1,\ldots n-1)(i=1,2\ldots n). \]

Formula (1-9) determines the approximate value of the coefficients in the expansion (1-4), obtained as a result of the interpolation process at the zeros of the first of the discarded “\(e\)”-polynomials in (1-4), i.e., at the zeros of \(T_n^*(t)\). If we normalize the vectors \(T_k\) \((k=0,1,\ldots,n-1)\), then we obtain an orthonormal basis of the subspace \(R_n^{\tau\tau}\), which may be regarded as the set of eigenvectors of a certain self-adjoint operator \(A^\tau\). The system of equations for the normalized vectors

\[ \overline{T}_k=\sqrt{\frac{2}{n}}\cdot T_k \]

has the form

\[ A^\tau\cdot \overline{T}_k=\lambda_k\cdot \overline{T}_k \quad (k=0,1,\ldots,n-1), \tag{1-10} \]

where \(\lambda_k\) \((k=0,1,\ldots,n-1)\) are the eigenvalues of the operator \(A^\tau\). Introduce the \(n\times n\) matrix \(U_n^\tau\), whose columns are the coordinates of the vectors \(\overline{T}_k\) \((k=0,1,\ldots,n-1)\). Obviously,

\[ U_n^\tau\cdot \widetilde{U}_n^\tau=E, \tag{1-11} \]

where \(E\) is the identity matrix, and the symbol \((\sim)\) denotes transposition. It can also be seen that

\[ \widetilde{U}_n^\tau=\sqrt{\frac{n}{2}}\cdot \widetilde{I}_n^{\tau\tau}, \tag{1-12} \]

and, consequently,

\[ I_n^{\tau\tau}\cdot \widetilde{I}_n^{\tau\tau}=\frac{2}{n}E. \tag{1-13} \]

Multiply equation (1-7) by \(I_n^{\tau\tau}\); then we obtain

\[ f^\tau=\frac{n}{2}I_n^{\tau\tau}b^\tau \tag{1-14} \]

This equation is equivalent to the vector equation (1-5), obtained as a result of interpolation; therefore the matrix \(\frac{n}{2}I_n^{\tau\tau}\) should naturally be called the interpolation matrix. The columns of this matrix are the vectors \(T_k\) \((k=0,1,\ldots,n-1)\), multiplied by \(\frac{2}{n}\).

Let us find the structure of the operator \(A^\tau\), whose eigenvectors form the orthonormal basis of the subspace under consideration. Introduce the diagonal matrix \(\lambda\)

\[ \lambda= \begin{bmatrix} \lambda_0 & 0 & 0 & \ldots & 0\\ 0 & \lambda_1 & 0 & \ldots & 0\\ \ldots & \ldots & \ldots & \ldots & \ldots\\ 0 & 0 & 0 & \ldots & \lambda_{n-1} \end{bmatrix} \tag{1-15} \]

and write the matrix equation

\[ A^\tau\cdot U_n^\tau=U_n^\tau\lambda. \tag{1-16} \]

It is equivalent to the entire set of equations of the form (1-10) [4]. Multi-

Multiplying (1-16) on the right by the transposed matrix \(\widetilde U_n^{\,T}\). Taking into account (1-11), we shall have the following representation for the matrix operator

\[ A^{T}=U_n^{\lambda}\widetilde U_n^{\,T}=\frac{n}{2} I_n^{TT,\lambda}\widetilde I_n^{TT}. \tag{1-17} \]

b) The subspace \(R_n^{SS}\)

An orthogonal basis generating the subspace \(R_n^{SS}\) can be obtained analogously to the preceding one.

Take the first \(n\) “\(e\)” functions of Chebyshev of the third kind \(S_1(t),\ldots,S_n(t)\) [5]. Let \(t_i^S\) \((i=1,2,\ldots,n)\) be the zeros of the function \(S_{n+1}(t)\), i.e., the zeros of the “\(e\)” Chebyshev polynomial of the second kind \(U_n^*(t)\), since

\[ S_{n+1}(t)=2\cdot e^{-\frac{at}{2}}\sqrt{1-e^{-at}}\cdot U_n^*(t). \]

Form \(n\) vectors by the formula

\[ S_k=\{S_k(t_1^S),\,S_k(t_2^S),\ldots,S_k(t_n^S)\}\quad (k=1,2,\ldots,n). \tag{1-18} \]

It can be shown [Appendix 1] that

\[ \begin{aligned} (S_k,S_m)&=0 \qquad k\ne m,\\ (S_k,S_k)&=\|S_k\|^2=\frac{n+1}{2}, \end{aligned} \left\} \tag{1-19} \right. \]

therefore the system of vectors \(S_1,S_2,\ldots,S_n\) may be taken as an orthogonal basis of a certain subspace, which we shall denote by \(R_n^{SS}\).

Consider the problem of interpolation by means of the functions \(S_k(t)\) \((k=1,2,\ldots,n)\). We impose on the function \(f(t)\), \((0\le t\le \infty)\), the additional condition

\[ f(0)=f(\infty)=0. \tag{1-20} \]

The problem consists in determining the coefficients of the partial sum

\[ f(t)\simeq \sum_{k=1}^{n}\beta_k^S S_k(t), \tag{1-21} \]

which has the vector analogue

\[ f^S=\sum_{k=1}^{n}\beta_k^S S_k. \tag{1-22} \]

The meaning of equality (1-22) is that the vector \(f^S\), belonging to the subspace \(R_n^{SS}\), is expanded in the orthogonal basis of this subspace; therefore

\[ \beta_k^S=\frac{(f^S,S_k)}{\|S_k\|^2} =\frac{2}{n+1}\sum_{i=1}^{n} f(t_i^S)S_k(t_i^S) \quad (k=1,2\ldots n), \tag{1-23} \]

where \(f(t_i^S)\) are the values of the function \(f(t)\) at the interpolation nodes, i.e., the coordinates of the vector \(f^S\). Analogously to the preceding case, we can introduce the interpolation matrix \(\dfrac{n+1}{2} I_n^{SS}\), whose columns are the coordinates of the corresponding vectors of the orthogonal basis, multiplied by \(\dfrac{2}{n+1}\).

With the help of this matrix, (1-22) can be written in the form (1-14)

\[ f^S=\frac{n+1}{2} I_n^{SS}\beta^S, \tag{1-24} \]

where \(\beta^S=(\beta_1^S,\beta_2^S,\ldots,\beta_n^S)\) is a vector of the \(n\)-dimensional Euclidean space \(R_n\). Thus, \(\dfrac{n+1}{2} I_n^{SS}\) is a matrix effecting an isomorphic correspondence between the space \(R_n\) and its subspace \(R_n^{SS}\). Taking into account that

\[ \frac{n+1}{2}\, I_n^{SS}\cdot \widetilde I_n^{SS}=E, \]

one can obtain the solution of equation (1-24) with respect to the vector \(\beta^S\)

\[ \beta^S=I_n^{SS}\cdot f^S, \tag{1-25} \]

which combines the whole set of scalar equations (1-23). We also note that the orthonormal basis \(S_1,S_2,\ldots,S_n\) of the subspace under consideration consists of eigenvectors of the matrix operator \(A^S\), whose structure is determined by the formula

\[ A^S=\frac{n+1}{2}\, I_n^{SS}\cdot \lambda \cdot \widetilde I_n^{SS}. \]

c) The subspace \(R_n^{ST}\)

Let the finite expansion (1-21) hold. As interpolation nodes we choose not the zeros of the “\(e\)” polynomial \(U_n^*(t)\), as in the preceding case, but the zeros of the “\(e\)” polynomial \(T_n^*(t)\). The vector expansion will take the form

\[ f^T=\sum_{k=1}^{n}\beta_k^T\cdot S_k^T, \tag{1-26} \]

where \(f^T\) is (1-6), and the vector \(S_k^T\) is defined by the formula

\[ S_k^T=\{S_k(t_1^T),S_k(t_2^T),\ldots,S_k(t_n^T)\}. \tag{1-27} \]

The system of vectors \(S_k^T\) \((k=1,2,\ldots,n)\) also has the orthogonality property [Appendix 1]

\[ \begin{aligned} (S_k^T,S_m^T)&=0 \quad k\ne m,\\ (S_k^T,S_k^T)&=\|S_k^T\|^2=\frac{n}{2} \end{aligned} \tag{1-28} \]

and, consequently, can serve as an orthogonal basis generating a certain subspace \(R_n^{ST}\). All the formulas obtained earlier retain their form. In particular, if one introduces the interpolation matrix \(\dfrac{n}{2} I_n^{ST}\), whose columns are the coordinates of the corresponding vectors of the orthogonal basis, multiplied by \(\dfrac{n}{2}\), then (1-26) can be written in the form

\[ f^T=\frac{n}{2}\, I_n^{ST}\cdot \beta^T, \tag{1-29} \]

where

\[ \beta^T=(\beta_1^T,\beta_2^T,\ldots,\beta_n^T). \tag{1-30} \]

We shall also have

\[ \beta_k^T=\frac{(f^T,S_k^T)}{\|S_k^T\|^2} =\frac{2}{n}\sum_{i=1}^{n} f(t_i^T)\cdot S_k(t_i^T) \quad (k=1,2\ldots n) \tag{1-31} \]

or, more compactly,

\[ \beta^T=\widetilde I_D^{ST}\cdot f^T. \tag{1-32} \]

The simultaneous solution of (1-29) and (1-14) gives

\[ \beta^{\mathrm T}=\frac{n}{2}\,\widetilde I_n^{\mathrm{ST}}\cdot I_n^{\mathrm{TT}}\cdot b^{\mathrm T}, \]

\[ b^{\mathrm T}=-\frac{n}{2}\,\widetilde I_n^{\mathrm{TT}}\cdot I_n^{\mathrm{ST}}\cdot \beta^{\mathrm T}. \tag{1-33} \]

These vector equalities relate the coefficients of the expansions (1-4) and (1-21), found by interpolation at the same nodes for one and the same function. The matrices \(\frac{n}{2}\widetilde I_n^{\mathrm{ST}}\cdot I_n^{\mathrm{TT}}\) and \(\frac{n}{2}\widetilde I_n^{\mathrm{TT}}\cdot I_n^{\mathrm{ST}}\) are matrices transforming the vector \(\beta^{\mathrm T}\) into the vector \(b^{\mathrm T}\), and conversely.

2. Relation to Numerical Integration

The aim of the interpolation processes under consideration is to obtain approximate values of the Fourier coefficients in finite expansions in terms of the “\(e\)”-functions \(T_k^*(t)\) and \(S_k(t)\) \((k=0,1,\ldots,n)\). The Fourier coefficients of these expansions are certain integrals [5]; therefore the interpolation formulas obtained are formulas of numerical quadrature. We shall show that they belong to the group of quadrature formulas of highest accuracy, i.e., to the so-called Gaussian quadratures [3].

Consider the integral

\[ \frac{2a}{\pi}\int_0^\infty \frac{\varphi(t)e^{-\frac{at}{2}}}{\sqrt{1-e^{-at}}}\,dt. \tag{2-1} \]

Make the substitution \(x=e^{-at}\); then we obtain

\[ \frac{2}{\pi}\int_0^1 \frac{\varphi^*(x)}{\sqrt{x(1-x)}}\,dx. \tag{2-2} \]

The weight function

\[ P(x)=\frac{1}{\sqrt{x(1-x)}} \]

is a special case of the Jacobi weight function

\[ P^{(\gamma\delta)}(x)=x^\gamma(1-x)^\delta, \]

corresponding to the values \(\gamma=\delta=-\frac{1}{2}\). This special case is characterized by the fact that the Jacobi polynomials are transformed into the shifted Chebyshev polynomials \(T_n^*(x)\) \((n=0,1,2,\ldots)\), which form a complete orthogonal system of functions on \((0,1)\) with weight \(P^{\left(-\frac12,-\frac12\right)}(x)\) [3, 4]. According to the general theory of Gaussian quadratures [3], an approximate value of the integral (2-2) can be obtained from \(n\) ordinates of the function \(\varphi^*(x)\) at the points \(x_i^{\mathrm T}\) \((i=1,2,\ldots,n)\), where

\[ x_i^{\mathrm T}=\cos^2\frac{(2i-1)\pi}{4n} \qquad (i=1,2,\ldots,n) \tag{2-3} \]

are the zeros of the polynomial \(T_n^*(x)\) [Appendix 1]. Thus,

\[ \frac{2}{\pi}\int_0^1 \frac{\varphi^*(x)}{\sqrt{x(1-x)}}\,dx = \frac{2}{n}\sum_{i=1}^n \varphi^*(x_i^{\mathrm T})+R_n(\varphi^*), \tag{2-4} \]

where \(R_n(\varphi^*)\) is the remainder term, having the form [3]

\[ R_n(\varphi^*)=\frac{1}{2^{2n-2}}\cdot \frac{\varphi^{*(2n)}(\xi)}{(2n)!},\qquad 0<\xi<1. \tag{2-5} \]

Let us return to the former variable; then, obviously,

\[ \frac{2a}{\pi}\int_{0}^{\infty} \frac{e^{-\frac{at}{2}}\varphi(t)}{\sqrt{1-e^{-at}}}\,dt = \frac{2}{n}\sum_{i=1}^{n}\varphi(t_i^T)+R_n^T(\varphi). \tag{2-6} \]

Here \(t_i^T\) \((i=1,2,\ldots,n)\) are the zeros of the Chebyshev “\(e\)-polynomial of the first kind \(T_n^*(t)\) [Appendix 1]. Put

\[ \varphi(t)=\varphi_T(t)=f(t)T_k^*(t)\qquad (k=0,1,2,\ldots,n-1), \tag{2-7} \]

then the integral (2-1) will determine the Fourier coefficients of the function \(f(t)\)

\[ b_k=\frac{2a}{\pi}\int_{0}^{\infty} \frac{e^{-\frac{at}{2}}f(t)T_k^*(t)}{\sqrt{1-e^{-at}}}\,dt \qquad (k=0,1\ldots n-1), \tag{2-8} \]

i.e., the coefficients in the finite expansion of the function \(f(t)\) in Chebyshev “\(e\)-polynomials of the first kind [5]

\[ f(t)=\frac{b_0}{2}+\sum_{k=1}^{n-1} b_k T_k^*(t), \tag{2-9} \]

According to (2-6), taking (2-7) and (2-8) into account, we obtain

\[ b_k=\frac{2a}{\pi}\int_{0}^{\infty} \frac{e^{-\frac{at}{2}}f(t)}{\sqrt{1-e^{-at}}}\cdot T_k^*(t)\,dt = \frac{2}{n}\sum_{i=1}^{n} f(t_i^T)T_k^*(t_i^T)+R_n^T(\varphi_T). \tag{2-10} \]

If one has in view (1-9), then one may write

\[ b_k=b_k^T+R_n^T(\varphi_T), \tag{2-11} \]

where

\[ b_k^T=\frac{2}{n}\sum_{i=1}^{n} f(t_i^T)T_k^*(t_i^T) \tag{1-9} \]

is an approximate value of the Fourier coefficient in the expansion (2-9), obtained as a result of interpolation at the zeros of the “\(e\)-polynomial \(T_n^*(t)\).

The remainder term \(R_n^T(\varphi_T)\) has the form

\[ R_n^T(\varphi_T)= \frac{1}{2^{2n-2}}\cdot \frac{\dfrac{d^{2n}}{dt^{2n}}\left[f(t)T_k^*(\tau)\right]}{(2n)!}, \qquad 0<\tau<\infty. \tag{2-12} \]

Put now

\[ \varphi(t)=\varphi_S(t)=f(t)\cdot S_k(t)\qquad (k=1,2,\ldots,n), \tag{2-13} \]

and we shall assume that \(f(0)=f(\infty)=0\). This requirement does not reduce generality, since for \(f(0)\ne0\) and \(f(\infty)\ne0\) one may consider the function

\[ f_1(t)=f(t)-f(0)e^{-\frac{at}{2}}-f(\infty)\left(1-e^{-\frac{at}{2}}\right). \]

The case \(f(0)=\infty\) and \(f(\infty)=\infty\) is excluded.

\(S_k(t)\) \((k=1,2,\ldots,n)\) are the Chebyshev “\(e\)-functions of the third kind [5]. Substituting (2-13) into (2-1), the values of the integrals

\[ \beta_k=\frac{2a}{\pi}\int_{0}^{\infty} \frac{e^{-\frac{at}{2}}f(t)}{\sqrt{1-e^{-at}}}\,S_k(t)\,dt \qquad (k=1,2\ldots n) \tag{2-14} \]

will give the Fourier coefficients of the expansion of the function \(f(t)\) in a finite series in the orthogonal functions \(S_k(t)\) [5]

\[ f(t)\cong \sum_{k=1}^{n}\beta_k S_k(t). \tag{1-21} \]

If we substitute (2-13) into the quadrature formula (2-6), we obtain

\[ \beta_k=\frac{2}{n}\sum_{i=1}^{n} f(t_i^T)\cdot S_k(t_i^T)+R_n^T(\varphi_s) \qquad (k=1,2,\ldots,n). \]

According to (1-31)

\[ \beta_k^T=\frac{2}{n}\sum_{i=1}^{n} f(t_i^T)\cdot S_k(t_i^T) \qquad (k=1,2,\ldots,n) \tag{1-31} \]

is an approximate value of the Fourier coefficients in the expansion (1-21), obtained as a result of interpolation at the zeros of the “\(e\)”-polynomial \(T_n^*(t)\).

The quadrature formula (2-6) makes it possible to obtain an approximate value of the integral from \(n\) ordinates of the function \(\varphi(t)\) taken at the zeros of the “\(e\)”-polynomial \(T_n^*(t)\). This result is a consequence of the orthogonality of the polynomials \(T_k^*(t)\) \((k=0,1,2,\ldots,n)\) with weight

\[ \frac{e^{-\frac{a}{2}t}}{\sqrt{1-e^{-at}}}=P(t), \]

but with the same weight there are also orthogonal “\(e\)”-Chebyshev functions of the third kind \(S_k(t)\) \((k=1,2,\ldots)\), which are likewise eigenfunctions of the corresponding Sturm–Liouville differential operator [1, 4]. This circumstance makes it possible, without changing anything in essence, to obtain a quadrature formula from the ordinates of the function \(\varphi(t)\) taken at the zeros of the function \(S_{n+1}(t)\), i.e., at the points \(t_i^S\) \((i=1,2,\ldots,n+1)\) [Appendix 1]. We shall have \(n+2\) ordinates \(\varphi(t_i^S)\) \((i=0,1,2,\ldots,n+1)\), including the values at the endpoints of the interval. If the natural requirement \(\varphi(0)=\varphi(\infty)=0\) is satisfied, then the quadrature formula takes the form

\[ \frac{2a}{\pi}\int_{0}^{\infty} \frac{e^{-\frac{a}{2}t}\varphi(t)}{\sqrt{1-e^{-at}}}\,dt = \frac{2}{n+1}\sum_{i=1}^{n}\varphi(t_i^S)+R_n^S(\varphi). \tag{2-15} \]

Let \(\varphi(t)\) be defined by formula (2-13); then on the left we obtain the exact value of the Fourier coefficient \(\beta_k\). Thus,

\[ \beta_k=\frac{2}{n+1}\sum_{i=1}^{n} f(t_i^S)\cdot S_k(t_i^S)+R_n^S(\varphi_S) \tag{2-16} \]

or, in another form,

\[ \beta_k=\beta_k^S+R_n^S(\varphi_S), \]

where \(\beta_k^S\) is an approximate value of the Fourier coefficient, determined by formula (1-23):

\[ \beta_k^S=\frac{2}{n+1}\sum_{i=1}^{n} f(t_i^S)\cdot S_k(t_i^S). \tag{1-23} \]

3. Errors and comparison of interpolation processes

The accuracy of the interpolation processes under consideration is determined by the absolute value of the remainder terms in the corresponding

quadrature formulas. However, it is practically impossible to compare accuracy on this basis, since the value of the \(2n\)-th derivative of the function being expanded is required at some unknown point located in the interval \((0,\infty)\).

Below we consider another approach to this problem, one very convenient for purposes of comparison.

Let \(f(t)\) \((0\leq t\leq \infty)\) be a continuous and bounded function on the whole interval; then its expansion in an infinite series in the Chebyshev “\(e\)” polynomials \(T_m^*(t)\) will converge uniformly and represent \(f(t)\) exactly at every point:

\[ f(t)=\frac{b_0}{2}+\sum_{m=1}^{\infty} b_m T_m^*(t), \tag{3-1} \]

The first \(n\) coefficients of this expansion can be determined approximately by interpolating at the zeros of the “\(e\)” polynomial \(T_n^*(t)\), i.e., by formula (1-9)

\[ b_k^T=\frac{2}{n}\sum_{i=1}^{n} f(t_i^T)\cdot T_k^*(t_i^T) \quad (k=0,1,\ldots,n-1). \tag{1-9} \]

Put \(t=t_i^T\) in (3-1) and substitute into (1-9). Then we obtain

\[ b_k^T=\frac{b_0}{n}\sum_{i=1}^{n} T_k^*(t_i^T) +\frac{2}{n}\sum_{m=1}^{\infty} b_m \sum_{i=1}^{n} T_m^*(t_i^T)\cdot T_k^*(t_i^T). \]

If we take into account [Appendix 1] that

\[ \frac{2}{n}\sum_{i=1}^{n}T_m^*(t_i^T)\cdot T_k^*(t_i^T) = \frac{1}{2n}\left[ \frac{\sin(m-k)\pi}{\sin\frac{(m-k)\pi}{2n}} + \frac{\sin(m+k)\pi}{\sin\frac{(m+k)\pi}{2n}} \right], \tag{3-2} \]

then we may write

\[ b_k^T= \frac{b_0}{2}\sum_{i=1}^{n}T_k^*(t_i^T) + \frac{1}{2n}\sum_{m=1}^{\infty} b_m \left[ \frac{\sin(m-k)\pi}{\sin\frac{(m-k)\pi}{2n}} + \frac{\sin(m+k)\pi}{\sin\frac{(m+k)\pi}{2n}} \right]. \tag{3-3} \]

Resolving the indeterminacies for \(m=\nu 2n\pm k\) \((\nu=0,1,2,\ldots)\), we shall have (for all other values of \(m\), (3-3) is equal to zero)

\[ b_k^T=b_k+\sum_{\nu=1}^{\infty}\left(b_{\nu 2n-k}+b_{\nu 2n+k}\right) \quad (k=0,1\ldots,n-1) \tag{3-4} \]

or, in expanded form,

\[ b_k^T=b_k+b_{2n-k}+b_{4n-k}+\ldots+b_{2n+k}+b_{4n+k}+\ldots . \tag{3-5} \]

The expression

\[ \sum_{\nu=1}^{\infty}\left(b_{\nu 2n-k}+b_{\nu 2n+k}\right). \tag{3-6} \]

is, obviously, the error of the interpolation under consideration. For small \(k\) (in comparison with \(n-1\)) the error is small; however, as \(k\) increases the error grows, and for values of \(k\) close to \(n-1\) it may be commensurate with the value of \(b_k^T\). Thus, for \(n=7\), we obtain the following picture:

\[ \begin{aligned} b_0^T&=b_0+2b_{14}+2b_{28}+\ldots,\\ b_1^T&=b_1+b_{13}+b_{15}+\ldots,\\ &\cdot\quad \cdot\quad \cdot\quad \cdot\quad \cdot\quad \cdot\\ b_5^T&=b_5+b_9+b_{19}+\ldots,\\ b_6^T&=b_6+b_8+b_{20}+\ldots . \end{aligned} \]

If the function \(f(t)\) is sufficiently smooth, i.e., the series (3-1) converges rapidly, then the interpolation error will be small. Let us impose an additional restriction on the function \(f(t)\) being expanded. We require that the equality

\[ f(0)=f(\infty)=0, \tag{3-7} \]

be satisfied; then the expansion in the orthogonal system \(S_m(t)\) \((m=1,2,\ldots)\) will converge uniformly:

\[ f(t)=\sum_{m=1}^{\infty}\beta_m S_m(t). \tag{3-8} \]

If interpolation at the zeros of the “\(e\)-polynomial” \(T_n^*(t)\) is used, then the approximate value of the first \(n\) coefficients of the expansion (3-8) can be determined by formula (1-31):

\[ \beta_k^T=\frac{2}{n}\sum_{i=1}^{n} f(t_i^T)S_k(t_i^T)\quad (k=1,2,\ldots,n). \tag{1-31} \]

Proceeding analogously to the preceding case, we obtain

\[ \beta_k^T=\frac{2}{n}\sum_{m=1}^{\infty} b_m \sum_{i=1}^{n} S_m(t_i^T)\cdot S_k(t_i^T). \]

Taking into account that [Appendix 1]

\[ \frac{2}{n}\sum_{i=1}^{n} S_m(t_i^T)S_k(t_i^T) = \frac{1}{2n} \left[ \frac{\sin(m-k)\pi}{\sin \dfrac{(m-k)\pi}{2n}} - \frac{\sin(m+k)\pi}{\sin \dfrac{(m+k)\pi}{2n}} \right], \]

we shall have

\[ \beta_k^T= \frac{1}{2n}\sum_{m=1}^{\infty}\beta_m \left[ \frac{\sin(m-k)\pi}{\sin \dfrac{(m-k)\pi}{2n}} - \frac{\sin(m+k)\pi}{\sin \dfrac{(m+k)\pi}{2n}} \right]. \tag{3-9} \]

All terms of this series for which \(m\ne \nu 2n\pm k\) \((\nu=0,1,\ldots)\) will be equal to zero. Resolving the indeterminacies for \(m=\nu 2n\pm k\) \((\nu=0,1,\ldots)\), we obtain

\[ \beta_k^T=\beta_k+\beta_{2n-k}-\beta_{4n-k}+\ldots-\beta_{2n+k}+\beta_{4n+k}-\ldots . \tag{3-10} \]

All the remarks made earlier remain valid for this case as well. There is, however, one peculiarity, namely that the last coefficient, corresponding to \(k=n\), can be determined as a result of interpolation quite accurately. Indeed, putting \(k=n\) in (3-10), we have

\[ \beta_n^T=2\beta_n-2\beta_{3n}+\ldots, \]

whence, with high accuracy, one may take

\[ \frac{1}{2}\beta_n^T=\beta_n. \]

Thus, if the last coefficient is taken with weight \(\dfrac{1}{2}\), then the value obtained will differ very little from the Fourier coefficient \(\beta_n\). In other words, the desired expansion should be written in the form

\[ f(t)\cong \beta_1^T S_1(t)+\beta_2^T S_2(t)+\ldots+\frac{1}{2}\beta_n^T S_n(t), \tag{3-11} \]

if the coefficients are determined by formula (1-31).

Let us find the relationship between the coefficient \(\beta_k^S\), obtained by interpolation at the zeros of the function \(S_{n+1}(t)\), and the Fourier coefficients in the expansion (3-8).

We have

\[ \beta_k^S=\frac{2}{n+1}\sum_{i=1}^n f(t_i^S)S_k(t_i^S)\qquad (k=1,2,\ldots,n). \tag{1-23} \]

Repeating the same procedure, we obtain [see Appendix 1]

\[ \beta_k^S=\frac{1}{2(n+1)}\sum_{m=1}^{\infty}\beta_m \left[ \frac{\sin(m-k)\pi\cdot \cos\dfrac{(m-k)\pi}{2(n+1)}}{\sin\dfrac{(m-k)\pi}{2(n+1)}}- \right. \]

\[ \left. -\frac{\sin(m+k)\pi\cdot \cos\dfrac{(m+k)\pi}{(n+1)}}{\sin\dfrac{(m+k)\pi}{2(n+1)}} \right], \tag{3-12} \]

whence, resolving the indeterminacies for \(m=\nu 2(n+1)\pm k\) \((\nu=0,1,\ldots)\), we find the desired relation

\[ \beta_k^S=\beta_k-\beta_{2(n+1)-k}-\beta_{4(n+1)-k}-\ldots +\beta_{2(n+1)+k}+\beta_{4(n+1)+k}+\ldots . \tag{3-13} \]

The maximum absolute error in such interpolation will be equal to

\[ |\beta_k^S-\beta_k|= |\beta_{2(n+1)-k}|+|\beta_{4(n+1)-k}|+\ldots+ |\beta_{2(n+1)+k}|+\ldots . \tag{3-14} \]

Whereas, in interpolation at the zeros of the “\(e\)”-polynomial \(T_n^*(t)\), it has the form

\[ |\beta_k^T-\beta_k|= |\beta_{2n-k}|+|\beta_{4n-k}|+\ldots+ |\beta_{2n+k}|+\ldots . \tag{3-15} \]

If the coefficients of the expansion (3-8), beginning with \(m>n\), decrease monotonically in absolute value, then

\[ |\beta_k^S-\beta_k|<|\beta_k^T-\beta_k|, \]

i.e., the coefficients determined by formula (1-23) differ less from the corresponding Fourier coefficients than do the coefficients \(\beta_k^T\), determined by the interpolation formula (1-31), for one and the same number of ordinates. In other words, interpolation at the zeros of \(S_{n+1}(t)\) is appreciably more accurate than interpolation at the zeros of \(T_n^*(t)\).

To obtain the same accuracy, the coefficients \(\beta_k^T\) must already be computed from \(n+1\) ordinates of the function \(f(t)\), since replacing \(n\) by \(n+1\) in (3-15) gives (3-14).

As for the coefficients \(b_k^T\), determined by the interpolation formula (1-9), their difference from the corresponding Fourier coefficients in the expansion (3-1) will be determined by a formula of the type (3-15); however, the last coefficient \(b_{n-1}^T\) will have a significant error. Therefore such interpolation will, generally speaking, be less accurate than when the functions \(S_k(t)\) are used.

The latter can be explained by the fact that, when these functions are used, the zero boundary conditions are satisfied in advance, independently of the expansion coefficients; therefore the equality in (3-11) is in fact attained not at \(n\) points of the interval \((0,\infty)\), but at \(n+2\) points, which cannot but entail some increase in accuracy.

Thus, of all the interpolation processes considered, the greatest accuracy is possessed by the process by which the coefficients \(\beta_k^S\) \((k=1,2,\ldots,n)\) are determined, i.e., the quadrature formula (2-15).

Let us also note that, from the uniform convergence of the expansions (3-1) and (3-8), it follows that \(b_m\to 0\) and \(\beta_m\to 0\) as \(m\to\infty\). Taking this into account, one can

it can be seen from (3-5), (3-10), and (3-13) that, as \(n \to \infty\), the approximate values of the Fourier coefficients tend to the exact ones, i.e., the interpolation processes converge.

Appendix 2 gives the values of the zeros of the “\(e\)”-polynomial \(T_n^*(t)\) and of the function \(S_{n+1}(t)\) up to \(n=8\). Because of lack of space, the interpolation matrices are not given. Their elements can be found with the aid of tables for trigonometric functions. Thus, the elements of the matrix

\[ \frac{n}{2} I_n^T \]

will be the quantities [Appendix 1]

\[ T_k^*(t_i^T)=\cos \kappa x_i^T =\cos \kappa \frac{(2i-1)\pi}{2n} \qquad (\kappa=1\ldots n-1)\quad (i=1,2\ldots n). \]

For the elements of the matrices

\[ \frac{n}{2} I_n^{ST} \quad \text{and} \quad \frac{n+1}{2} I_n^{SS}, \]

respectively, we shall have the formulas

\[ S_k(t_i^T)=\sin \kappa \alpha_i^T =\sin \kappa \frac{(2i-1)\pi}{2n}, \]

\[ S_k(t_i^S)=\sin k \alpha_i^S =\sin k \frac{i\pi}{n+1}. \]

The columns of these matrices are the coordinates of orthogonal vectors forming bases of the corresponding subspaces. As for the parameter “\(a\),” it must be chosen from the conditions indicated in [6]. It is necessary, however, to remember that if \(f(0)\ne 0\) and \(f(\infty)\ne 0\), then the expansion coefficients with respect to the functions \(S_k(t)\) \((k=1,2,\ldots,n)\) should be sought for a new function of the form

\[ f_1(t)=f(t)-f(0)e^{-\frac{at}{2}}-f(\infty)\left(1-e^{-\frac{at}{2}}\right), \]

then for the given function we obtain the representation

\[ f(t)\simeq f(0)e^{-\frac{at}{2}} + f(\infty)\left(1-e^{-\frac{at}{2}}\right) + \sum_{k=1}^{n}\beta_k S_k(t). \]

As an example, let us find the approximate values of the Fourier coefficients \(b_k^T\), \(\beta_k^T\), and \(\beta_k^S\) for the function

\[ f(t)=e^{-t}\cdot \cos 3t. \tag{3-16} \]

This function has a rather strongly pronounced oscillatory character; however, we shall try to find the approximate values of several of the first coefficients, taking only eight ordinates in all. The results of the calculation by the schemes considered are given in Table 1 \((a=1)\).

Table 1.

The calculation was carried out using Chambers’ six-place mathematical tables; however, only four digits were taken, and in some cases rounding was performed. The last row gives the values of the Fourier coefficients \(\beta_k\) \((k=1,2,\ldots,6)\), obtained from the exact formulas [6]. Despite such a small number of ordinates for a function such as (3-16), the error is quite acceptable for engineering applications.

Appendix 1

Proof of the Orthogonality of the Basis Vectors

The substitution \(e^{-at}=\cos^2 \dfrac{\alpha}{2}\), or \(t=-\dfrac{2}{a}\ln\cos\dfrac{\alpha}{2}\) (P1-1), transforms the sequences of “e”-Chebyshev polynomials of the first kind \(T_k^*(t)\) \((k=0,1,2,\ldots)\) and “e”-functions \(S_k(t)\) \((k=1,2,\ldots)\) into sequences of trigonometric functions \(\cos k\alpha\) and \(\sin k\alpha\) \((k=0,1,2,\ldots)\), respectively [5].

The roots of the equation \(\cos n\alpha=0\) \((0\le \alpha\le \pi)\) are determined by the formula

\[ \alpha_i^T=(2i-1)\frac{\pi}{2n}\quad (i=1,2\ldots n), \tag{P1-2} \]

therefore the zeros of the “e”-polynomial \(T_n^*(t)\) are located at the points

\[ t_i^T=-\frac{2}{a}\ln\cos(2i-1)\frac{\pi}{4n}\quad (i=1,2\ldots n). \tag{P1-3} \]

For the equation \(\sin(n+1)\alpha=0\), on the open interval \((0,\pi)\) we shall have \(n\) roots:

\[ \alpha_i^S=\frac{i\pi}{n+1}\quad (i=1,2\ldots n). \tag{P1-4} \]

Accordingly, the zeros of the “e”-Chebyshev function \(S_{n+1}(t)\), coinciding with the zeros of the “e”-polynomial \(U_n^*(t)\), are located at the points

\[ t_i^S=-\frac{2}{a}\ln\cos\frac{i\pi}{2(n+1)}\quad (i=1,2\ldots n). \tag{P1-5} \]

a) Let us show that

\[ (T_k,T_m)=\sum_{i=1}^{n}T_k^*(t_i^T)T_m^*(t_i^T)= \begin{cases} 0, & m\ne k,\\[2mm] \dfrac{n}{2}, & m=k. \end{cases} \tag{P1-6} \]

Since

\[ T_k^*(t_i^T)=\cos k\alpha_i^T,\quad \text{and}\quad T_m^*(t_i^T)=\cos m\alpha_i^T, \]

the problem reduces to proving the relation

\[ \sum_{i=1}^{n}\cos k\alpha_i^T\cdot \cos m\alpha_i^T = \begin{cases} 0, & m\ne k,\\[2mm] \dfrac{n}{2}, & m=k. \end{cases} \]

We have

\[ \cos k\alpha_i^T\cdot \cos m\alpha_i^T = \frac{1}{2}\left[\cos(k-m)\alpha_i^T+\cos(k+m)\alpha_i^T\right]= \]

\[ = -\frac{1}{2}\cos(2i-1)\frac{(k-m)\pi}{2n} +\frac{1}{2}\cos(2i-1)\frac{(k+m)\pi}{2n}. \]

If we take into account that [2]

\[ \sum_{i=1}^{n}\cos(2i-1)x = \frac{1}{2}\frac{\sin nx}{\sin x}, \]

then we can write

\[ \sum_{i=1}^{n}\cos k\alpha_i^T \cdot \cos m\alpha_i^T = \frac{1}{4} \left[ \frac{\sin(k-m)\pi}{\sin \dfrac{(k-m)\pi}{2n}} + \frac{\sin(k+m)\pi}{\sin \dfrac{(k+m)\pi}{2n}} \right]. \]

For any unequal \(k\) and \(m\) this expression is equal to zero \((k+m\ne 2n\) and \(k-m\ne 2n)\). For \(k=m\), resolving the indeterminacy, we obtain \(\dfrac{n}{2}\). Thus, relation (P1-6) has been proved.

b) Let us show that

\[ (S_k,S_m)=\sum_{i=1}^{n}S_k(t_i^S)S_m(t_i^S) = \begin{cases} 0, & k\ne m,\\[4pt] \dfrac{n+1}{2}, & k=m. \end{cases} \tag{P1-} \]

Analogously to the preceding, the problem reduces to proving the relation

\[ \sum_{i=1}^{n}\sin k\alpha_i^S \cdot \sin m\alpha_i^S = \begin{cases} 0, & k\ne m,\\[4pt] \dfrac{n+1}{2}, & k=m. \end{cases} \]

Using (P1-4) and carrying out the transformation, we obtain

\[ \sum_{i=1}^{n}\sin k\alpha_i^S \cdot \sin m\alpha_i^S = \frac{1}{2}\sum_{i=1}^{n}\cos i\frac{(k-m)\pi}{n+1} - \frac{1}{2}\sum_{i=1}^{n}\cos i\frac{(k+m)\pi}{n+1}. \]

For any \(x\) the identity [2] holds

\[ \sum_{i=1}^{n}\cos ix = \frac{\cos \dfrac{n+1}{2}x \cdot \sin n\dfrac{x}{2}}{\sin \dfrac{x}{2}} = \frac{1}{2}\frac{\sin 2(n+1)x \cdot \cos \dfrac{x}{2}}{\sin \dfrac{x}{2}} - \cos^2 \frac{n+1}{2}x. \]

Substituting \(x \to \dfrac{(k-m)\pi}{n+1}\) and \(x \to \dfrac{(k+m)\pi}{n+1}\), we obtain, after simple transformations,

\[ \sum_{i=1}^{n}\sin k\alpha_i^S \cdot \sin m\alpha_i^S = \frac{1}{2} \left[ \frac{ \dfrac{1}{2}\sin(k-m)\pi \cdot \cos \dfrac{(k-m)\pi}{2(n+1)} }{ \sin \dfrac{(k-m)\pi}{2(n+1)} } - \frac{ \dfrac{1}{2}\sin(k+m)\pi \cdot \cos \dfrac{(k+m)\pi}{2(n+1)} }{ \sin \dfrac{(k+m)\pi}{2(n+1)} } - \sin k\pi \cdot \sin m\pi \right]. \]

For \(k\ne m\) this expression is equal to zero; for \(k=m\), resolving the indeterminacy, we obtain \(\dfrac{n+1}{2}\), which proves (III-7).

c) Let us show that

\[ (S_k^\tau, S_m^\tau)=\sum_{i=1}^{n} S_k(t_i^\tau)\cdot S_m(t_i^\tau) = \begin{cases} 0, & k\ne m,\\[4pt] \dfrac{n}{2}, & k=m. \end{cases} \tag{III-8} \]

Substitution of (III-3) gives

\[ (S_k^\tau, S_m^\tau)= \sum_{i=1}^{n}\sin k\alpha_i^\tau\cdot \sin m\alpha_i^\tau = \frac{1}{4} \left[ \frac{\sin (k-m)\pi}{\sin \dfrac{(k-m)\pi}{2n}} - \frac{\sin (k+m)\pi}{\sin \dfrac{(k+m)\pi}{2n}} \right], \]

whence (III-8) follows.

Appendix 2.

Zeros of the polynomial \(T_n^*(t)\)

Zeros of \(S_{n+1}(t)\), \(0<t<\infty\)

REFERENCES

1. B. Z. Vulikh. Introduction to Functional Analysis. Fizmatgiz, Moscow, 1958.

2. I. S. Gradshtein, I. M. Ryzhik. Tables of Integrals, Sums, Series, and Products. Fizmatgiz, Moscow, 1962.

3. V. I. Krylov, L. T. Shulgina. Handbook of Numerical Integration. Nauka, Moscow, 1966.

4. K. Lanczos. Practical Methods of Applied Analysis. Fizmatgiz, Moscow, 1961.

5. V. M. Osipov. “Exponential Polynomials and the Expansion of Certain Standard Signals.” Proceedings of TPI, vol. 180, Tomsk, 1969.

6. V. M. Osipov. “On the Question of the Approximate Inversion of the Laplace Transform.” (This volume).