INTERPOLATION METHOD FOR THE ANALYSIS OF LINEAR TIME-VARYING SYSTEMS

V. M. Osipov, A. N. Barkovsky

(Presented by the scientific-technical seminar of the AVTF TPI)

On the basis of the interpolation method for solving differential equations, the problem of analyzing linear systems with variable parameters is considered.

The input and output signals of the system are represented in the form of “point” vectors, whose coordinates are the values of real continuous processes at certain points (interpolation nodes). The concept is introduced of a system transformation matrix establishing a simple algebraic relation between the point vectors of the input and output signals. It is shown that the basic questions in the analysis of time-varying systems can be solved by means of rather simple and uniform algorithms, convenient both when using digital computers and in hand calculations.

Interpolation method for solving differential equations

Consider the differential equation

\[ L y = \frac{d^n y}{d t^n} + c_1(t)\frac{d^{\,n-1}y}{d t^{\,n-1}} + \ldots + c_{n-1}(t)\frac{dy}{dt} + c_n(t)y = f(t) \tag{1} \]

with initial conditions

\[ \left.\frac{d^\nu y}{d t^\nu}\right|_{t=0} = y_0^{(\nu)};\quad (\nu = 0,1,2,\ldots,n-1). \tag{2} \]

An approximate solution of equation (1) is sought in the form

\[ y_n(t) = \psi_0(t) + \sum_{k=0}^{m} a_k \varphi_k(t), \tag{3} \]

where \(\psi_0(t)\) is a function satisfying the given initial conditions;
\(\varphi_0(t), \varphi_1(t), \ldots, \varphi_m(t)\) is a prescribed system of functions satisfying zero initial conditions.

The general idea of the interpolation method is that the coefficients \(a_k\) \((k=0,1,2,\ldots,m)\), which enter into the approximate solution (3), are found from the conditions that the residual \(L y_n - f\) vanish at certain points, called interpolation nodes. This leads to a system of linear algebraic equations.

\[ \left. \frac{d^n}{dt^n}\left[\psi_0(t)+\sum_{k=0}^{m} a_k\varphi_k(t)\right] \right|_{t=t_i} + \left\{ \sum_{\nu=1}^{n} c_\nu(t)\, \frac{d^{\,n-\nu}}{dt^{\,n-\nu}} \left[\psi_0(t)+\sum_{k=0}^{m} a_k\varphi_k(t)\right] \right\}_{t=t_i} = f(t_i); \]
\[ (i=1,2,3,\ldots,m+1). \tag{4} \]

It is natural to expect that the convergence and accuracy of the method will depend essentially on the choice of the system of functions in the approximate solution (3) and on the choice of the interpolation nodes. Proceeding from these considerations, as the coordinate system of functions \(\varphi_k(t)\) \((k=0,1,\ldots,m)\) we choose the system of orthogonal exponential Chebyshev polynomials [1], and as the interpolation nodes—the zeros of the first of the discarded polynomials. In this case the approximate solution (3) takes the form

\[ y_n(t)=Ae^{-at}+B(e^{-at}-1)+C(e^{-at}-1)^2+\cdots+ \]
\[ +E(e^{-at}-1)^{n-1}+(e^{-at}-1)^n\sum_{k=0}^{m} a_kT_k^*(t), \tag{5} \]

where \(A,B,C,\ldots,E\) are constants determined from the initial conditions;
\(T_k^*(t)\) are exponential Chebyshev polynomials of the first kind [1].

In what follows, for simplicity, we shall assume that the initial conditions (2) are zero. Then

\[ y_n(t)=(e^{-at}-1)^n\sum_{k=0}^{m} a_kT_k^*(t). \tag{6} \]

Substituting (6) into (1) and carrying out simple transformations, the resulting differential equation can be reduced to the form

\[ \varphi_0^{1}(t)\left[\sum_{k=0}^{m} a_kT_k^*(t)\right]^{(n)} + \varphi_1^{1}(t)\left[\sum_{k=0}^{m} a_kT_k^*(t)\right]^{(n-1)} +\cdots+ \]
\[ +\varphi_n^{1}(t)\left[\sum_{k=0}^{m} a_kT_k^*(t)\right] =f(t). \tag{7} \]

Next, following the general idea of the interpolation method, we must construct a system of linear algebraic equations for determining the unknown parameters \(a_0,a_1,\ldots,a_m\).

Using the properties of the exponential polynomials and taking into account their orthogonality in the sense of summation over the interpolation nodes [2], it can be shown that the desired system of algebraic equations in vector–matrix form will have the form

\[ (-a)^n M_0 I^{T}_{m+1}D^n A + (-a)^{n-1}M_1 I^{T}_{m+1}D^{n-1}A+\cdots+ \]
\[ +(-a)M_{n-1}I^{T}_{m+1}DA + M_n I^{T}_{m+1}A =F, \tag{8} \]

where

\[ M_0,\ldots,M_n \]

are diagonal matrices of the form

\[ M_k=\operatorname{diag}\{\varphi_k^{1}(t_1),\ldots,\varphi_k^{1}(t_{m+1})\}. \]

\(I^{T}_{m+1}\) is the so-called interpolation matrix of order \(m+1\), whose elements are determined by the expression

\[ i_{\mu\nu}=T_{\nu-1}^{*}(t_\mu); \qquad (\mu=1,2,\ldots,m+1)(\nu=1,2,\ldots,m+1), \]

\(D\) is a differentiation matrix of the form

\[ D= \begin{bmatrix} 0 & \sqrt{2} & 2\sqrt{2} & \cdot & \cdot & \cdot & \cdot & m\sqrt{2}\\ 0 & 1 & 2\cdot 2 & \cdot & \cdot & \cdot & \cdot & 2m\\ 0 & 0 & 2 & \cdot & \cdot & \cdot & \cdot & 2m\\ \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot\\ 0 & \cdot & \cdot & \cdot & \cdot & \cdot & m-1 & 2m\\ 0 & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & m \end{bmatrix}. \]

\(A\) is the column vector of the unknown coefficients

\[ A=\{a_0,\ a_1,\ \ldots,\ a_m\}, \]

\[ F=\{f(t_1),\ f(t_2),\ \ldots,\ f(t_{m+1})\}. \]

\(a\) is a parameter of the Chebyshev polynomials [1].

From (8) we easily find the unknown vector

\[ A=\left[(-a)^n M_0 I^{\mathrm{TT}}_{m+1}D^n+(-a)^{n-1}M_1 I^{\mathrm{TT}}_{m+1}D^{n-1}+\ldots+M_n I^{\mathrm{TT}}_{m+1}\right]^{-1}F . \tag{9} \]

Using (6) further, one can obtain an analytic expression for the approximate solution.

As direct calculations have shown, the inverse matrix of expression (9) proves to be unstable. In order to overcome this kind of ill-posedness, it has proved effective to use an overdetermined system of algebraic equations with subsequent processing by the least-squares method. In this case the vector-matrix equation will have the form

\[ \left[(-a)^n M'_0 I^{\mathrm{TT}}_{m+j+1,m+1}D^n +(-a)^{n-1}M'_1 I^{\mathrm{TT}}_{m+j+1,m+1}D^{n-1} +\ldots+ M'_n I^{\mathrm{TT}}_{m+j+1,m+1}\right]A=F', \]

where \(M'_0, M'_1,\ldots,M'_n\) are diagonal matrices of the form

\[ M'_k=\operatorname{diag}\{\varphi_k^1(t_1),\ \varphi_k^1(t_2),\ \ldots,\ \varphi_k^1(t_{m+j+1})\} \tag{10} \]

\[ F'=\{f(t_1),\ f(t_2),\ \ldots,\ f(t_{m+j+1})\}. \]

\(I^{\mathrm{TT}}_{m+j+1,m+1}\) is an interpolation matrix of order \(m+j+1\), in which the last \(j\) columns have been discarded (it is assumed that the overdetermination has been carried out for \(j\) equations).

If we denote

\[ \left[(-a)^n M'_0 I^{\mathrm{TT}}_{m+j+1,m+1}D^n +(-a)^{n-1}M'_1 I^{\mathrm{TT}}_{m+j+1,m+1}D^{n-1} +\ldots+ M'_n I^{\mathrm{TT}}_{m+j+1,m+1}\right]=B, \]

then the problem of solving equation (10) reduces to computing the vector \(A\) by the formula

\[ A=(\widetilde{B}B)^{-1}\widetilde{B}F', \tag{11} \]

where \(\widetilde{B}\) is the matrix equal to the transpose of the matrix \(B\).

System transformation matrix. Determination of the transformation matrix for various connections of elements

We shall assume that the behavior of the nonstationary dynamic system under investigation is described by a differential equation of the form (1). Then the problem of determining the system response to the input signal \(f(t)\) reduces to solving this equation. In the interpolation method of solution, we represent the input signal of the system in the form of a vector, whose coordina-

components of which are the values of the function \(f(t)\) at the interpolation nodes (see (10)),

\[ F'=\{f(t_1), f(t_2), \ldots, f(t_{m+j+1})\}. \tag{12} \]

Using the interpolation method [2], from the found vector of coefficients (11) of the solution (6) it is not difficult to obtain the vector of approximate values of the system output signal at the same interpolation nodes, where

\[ Y_n=\eta I^{TT}_{m+j+1,m+1}A =\eta I^{TT}_{m+j+1,m+1}(\tilde{B}B)^{-1}\tilde{B}F', \tag{13} \]

\(\eta\) is a diagonal matrix which, in accordance with (6), has the form

\[ \eta=\operatorname{diag}\{(e^{-at_1}-1)^n,\ (e^{-at_2}-1)^n,\ \ldots,\ (e^{-at_{m+j+1}}-1)^n\}, \]

\(Y_n\) is the vector of the system output signal

\[ Y_n=\{y_n(t_1), y_n(t_2), \ldots, y_n(t_{m+j+1})\}. \tag{14} \]

Denoting

\[ \eta I^{TT}_{m+j+1,m+1}(\tilde{B}B)^{-1}\tilde{B}=W, \]

we finally find

\[ Y_n=WF'. \tag{15} \]

We shall call the matrix \(W\) the transformation matrix of the system.

Thus, if for a given dynamic system (1) the input signal is taken to be the vector (12), and the output signal the vector (14), then, according to (15), there exists a matrix \(W\) establishing the relation between the input and output signals of the system. This matrix, as is evident, does not depend on the applied excitation, but is completely determined by the structure of the original differential equation describing the system and, consequently, characterizes the properties of the system itself. Knowing the transformation matrix, one can, by means of the elementary operation of matrix multiplication, find the vector of approximate values of the response of the system under study to any input signal specified in the form of the point vector (12).

The transformation matrices of various connections of elements can be found from the transformation matrices of the elements by means of simple algebraic operations. Thus, for the series connection of \(m\) elements having transformation matrices \(W_1, W_2, \ldots, W_m\),

\[ W=W_m W_{m-1}\ldots W_1, \tag{16} \]

and, because of the noncommutativity of matrix multiplication, permutation of the matrices in the resulting expression is inadmissible.

For the parallel connection of \(m\) elements with transformation matrices \(W_1, W_2, \ldots, W_m\),

\[ W=W_1+W_2+\ldots+W_m. \tag{17} \]

The transformation matrix of a system encompassed by feedback is determined by the formula

\[ W=(E+W_1\cdot W_2)^{-1}W_1, \tag{18} \]

where

\[ \begin{aligned} W_1 &\text{ is the transformation matrix of the forward path;}\\ W_2 &\text{ is the transformation matrix of the feedback element;}\\ E &\text{ is the identity matrix.} \end{aligned} \]

Expressions (16), (17), (18) are reminiscent of analogous relations for systems with constant parameters.

Determination of the Probabilistic Characteristics of the Output Processes of Nonstationary Systems

As follows from what was set forth above, for a known transformation matrix of the system described by equation (1), the output random vector \(Y_n\) can be found from the input random vector \(F'\) by formula (15).

To determine the correlation function of the output signal of the system, we first reconstruct it as a function of time from its point vector, using the interpolation process.

The coefficients of the Fourier–Chebyshev series are found from the expression [2]

\[ A=\frac{2}{m+j+1}\,\widetilde{I}^{\mathrm{T}}_{m+j+1}Y_n =\frac{2}{m+j+1}\,\widetilde{I}^{\mathrm{T}}_{m+j+1}WF'. \tag{19} \]

Then the random function of the output process of the system can be represented approximately in the form

\[ y(t)=\widetilde{T}(t)A =\widetilde{T}(t)\frac{2}{m+j+1}\,\widetilde{I}^{\mathrm{T}}_{m+j+1}WF'_t, \tag{20} \]

where \(\widetilde{T}(t)\) is a row vector of Chebyshev polynomials

\[ \widetilde{T}(t)=\{T^*_0(t),T^*_1(t),\ldots,T^*_{m+j}(t)\}. \]

The index \(t\) on the vector \(F'\) means that the random function \(f(t)\), from which the vector \(F'\) is formed, is also a function of the argument \(t\). Similarly we find

\[ y(\tau)=\widetilde{T}(\tau)\frac{2}{m+j+1}\,\widetilde{I}^{\mathrm{T}}_{m+j+1}WF'_\tau \tag{21} \]

or, equivalently,

\[ y(\tau)=\frac{2}{m+j+1}\,\widetilde{F}'_\tau \widetilde{W}\,\widetilde{I}^{\mathrm{T}}_{m+j+1}T(\tau), \tag{22} \]

where

\[ I(\tau) \]

is a column vector of Chebyshev polynomials,

\[ \widetilde{F}'(\tau) \]

is a row vector.

Now form the product

\[ y(t)y(\tau)= \widetilde{T}(t)\frac{2}{m+j+1}\,\widetilde{I}^{\mathrm{T}}_{m+j+1}WF'_t \frac{2}{m+j+1}\,\widetilde{F}'_\tau\widetilde{W}\widetilde{I}^{\mathrm{T}}_{m+j+1}T(\tau) \]

and average it over the ensemble of realizations:

\[ \overline{y(t)y(\tau)} = \frac{4}{(m+j+1)^2}\, \widetilde{T}(t)\widetilde{I}^{\mathrm{T}}_{m+j+1}W\cdot \overline{F'_t\widetilde{F}'_\tau}\, \widetilde{W}\widetilde{I}^{\mathrm{T}}_{m+j+1}T(\tau). \tag{23} \]

The left-hand side of the resulting expression is the correlation function of the output process of the system under study, \(K_y(t,\tau)\). The averaged product of the column vector \(F'_t\) by the row vector \(\widetilde{F}'_\tau\) is a matrix of correlation moments of the form

\[ K_f(t_i,\tau_i)= \left| \begin{array}{cccc} K_f(t_1,\tau_1) & K_f(t_1,\tau_2) & \cdots & K_f(t_1,\tau_{m+j+1})\\ K_f(t_2,\tau_1) & K_f(t_2,\tau_2) & \cdots & K_f(t_2,\tau_{m+j+1})\\ \cdots & \cdots & \cdots & \cdots\\ K_f(t_{m+j+1},\tau_1) & K_f(t_{m+j+1},\tau_2) & \cdots & K_f(t_{m+j+1},\tau_{m+j+1}) \end{array} \right| \]

Taking this into account, we write expression (23) in the following form:

\[ K_y(t,\tau)=\frac{4}{(m+j+1)^2}\,\tilde{T}(t)\bar{I}^{\mathrm{tt}}_{m+j+1} W K_f(t_i,\tau_i)\bar{W}\bar{I}^{\mathrm{tt}}_{m+j+1}T(\tau). \tag{24} \]

To find the mathematical expectation of the output random process, it is sufficient to integrate, by the interpolation method, the original differential equation of system (1) with \(f(t)=m_f(t)\)—the mathematical expectation of the input random signal.

It can be shown that determination of the statistical characteristics of the system output process caused by the application of several random actions is performed by formulas analogous to (24).

REFERENCES

1. V. M. Osipov. Exponential polynomials and the expansion of certain typical signals. Izv. TPI, vol. 180. Tomsk, 1970.

2. V. M. Obshlov. Approximation of functions of time by the interpolation method. Izv. TPI, vol. 191. Tomsk, 1969.