A homogeneous high precision direct integration based on Chebyshev interpolation

1. Introduction

For the initial value problem of homogeneous linear differential equations with constant coefficients, the matrix index precise calculation method proposed by Academician Zhong Wanxie [1] can increase the calculation of matrix index to computer precision by selecting the integral step size within a reasonable range. This method is widely used in many fields such as random vibration, dynamic response and optimal control. For the processing of non-homogeneous term, Lin Jiahao et al. [2] obtained the analytical solution of Duhamel integral term for the non-homogeneous term of polynomial, periodic function, exponential function and their combination. However, the inverse matrix of the coefficient matrix of the equations is required to be solved. When the coefficient matrix is singular or nearly singular, this method cannot be used, thus greatly limiting its application range and reliability. In order to avoid the inverse operation of the coefficient matrix, Wang Mengfu et al. [3] directly calculates the Duhamel integral by numerical methods. Although such methods can achieve higher precision, the matrix index on the integration point needs to be additionally calculated. Xiang Yu et al. [4] solved the problem by polynomial approximation or Fourier series expansion of inhomogeneous terms. Such methods cannot carry out precision design. If the number of fitting times is higher, the Runge phenomenon will occur. Moreover, as the step size of the integration increases, a large error is generated, so the step size needs to be divided finely, and the amount of calculation is increased.

Based on the Chebyshev’s interpolation theory, this paper interpolates the sampling points of non-homogeneous terms. Different from the traditional uniform sampling, this paper can improve the calculation accuracy based on the Chebyshev interpolation point and avoid the Runge phenomenon. Using the error definition formula, the fitting times can be selected according to different precision requirements, which has strong controllability. Moreover, this method is easy to program and calculate on a computer, and the implementation is simple.

2. Precise calculation of matrix index

Differential equations in vibration engineering, structural dynamic response, etc. can all be transformed into state vector forms:

where $v\left(t\right)$ is the introduced state vector, $H$ is the constant matrix, and $r\left(t\right)$ is the structural load or control vector. Considering that when $r\left(t\right)=$ 0, the Eq. (1) is transformed into a homogeneous equation, and its solution is:

Let the integration step size $t-t_0=\tau$, then:

where $T=\mathrm{e}\mathrm{x}\mathrm{p}(H\bullet \tau )$. It can be known that the relationship between any two adjacent state vectors can be expressed as:

Therefore, the precise calculation of the transfer matrix becomes the core of the solution. In literature [1], the matrix index is processed as follows by using the addition theorem:

Let $m=2^N$ and the original integration step size become $∆t=\tau /m$. According to the literature [5], if $N=$ 20,$\tau =$ 0.01, the relative error is 0.151×10^-17, which has reached the full accuracy of the computer. This operation ensures the accuracy of the matrix index expansion with Taylor series. Take the first three items of expansion as follows:

where $T_a=H\bullet ∆t+{\left(H\bullet ∆t\right)}^2/\text{2}$, when calculating with a computer, since the value of $T_a$ is very small, it will become a mantissa when adding with $I$, and thus its original precision is lost. Therefore, it cannot be solved by direct power operation of $I+T_a$:

where $T_b=2T_a+T_a\times T_a$. So, calculating $T$ is equivalent to executing a loop statement:

The principle is shown in Fig. 1. Then execute the following statement at the end of the loop:

Fig. 1Schematic diagram of the loop

3. Chebyshev interpolation approximations

Polynomial approximation is applied to the non-homogeneous term of the differential equations, and the second-order polynomial is taken as an example:

when performing interpolation point sampling, if uniform sampling is performed as described in literature [4], the interpolation polynomial is liable to be twisted near the end of the interpolation interval. As the order increases, the twisting phenomenon becomes more pronounced and the error is greater. It’s the so-called Runge phenomenon.

The Chebyshev interpolation point is a specific optimal point spacing selection method that minimizes interpolation errors. Assuming that $n$ Chebyshev interpolation nodes are selected on the interval $[a$, $b]$, the node coordinates are respectively:

Then use the Newton’s divided difference formula [6] to find the values of the three vectors $r_0$, $r_1$, $r_2$. Extending the approximate polynomial to the original state vector of the differential equation:

Eq. (12) can be solved by the high precision direct integration method of the homogeneous equation described in Eq. (2) to Eq. (9).

4. Precision design

For Chebyshev interpolation, $f\left(x\right)$ is assumed to be the original function, $P\left(x\right)$ is a $n-1$ or lower order interpolation polynomial, then $n$ points $\left(x_1,y_1\right),\dots ,(x_n,y_n)$ are fitted, and the interpolation error is:

There is an inequality according to Chebyshev’s theory:

where the definition of $a$, $b$, $n$ is the same as Eq. (11). Then the maximum interpolation error can be written as:

For example, to perform a polynomial fitting on $\mathrm{s}\mathrm{i}\mathrm{n}x$ on [0, 1], it is required to be accurate to 6 decimal places. Through simple experiments, it is found that when $n=$ 7, the error bound ≈ 0.5715×10^-6, then the sixth-order polynomial fitting can be selected to meet the accuracy requirement.

5. Numerical examples

Differential equations with periodic function as non-homogeneous term:

Analytical solutions when non-homogeneous terms are sine/cosine functions were given in literature [2]:

where $r_1$ and $r_2$ are the coefficient vectors of the sine and cosine functions, respectively. Eq. (17) was verified by Mathematica and the results were consistent. Therefore, the result of Eq. (17) is taken as the exact solution. Introducing the state vector $V=\left\{x_1\left(t\right)x_2\left(t\right){\dot{x}}_1\left(t\right){\dot{x}}_2\right(t\left)\right\}$ and transforming the Eq. (16) into a first-order differential equation:

Then use the methods in literature [2], literature [1], literature [4] and this paper, respectively, to solve Eq. (18). The integration step $\tau =$ 0.01, the fine parameter $N=$ 20, the results of each method are shown in Table 1 and 2.

It can be seen from the Fig. 2 and Fig. 3 that because the literature [1] uses a linear approximation, so a large error will occur as the integration step size increases. The error of the literature [4] is close to the error of this paper. This is because the integration step size selected in the above calculation is small, and the larger integration step size is used to compare separately. Take $\tau =$ 0.015, repeat the above calculation process, and obtain the error curve comparison between the literature [4] and this paper are shown in Figs. 4, 5.

Fig. 2Error of x1

Fig. 3Error of x2

Fig. 4Error of x1 (τ= 0.015)

Fig. 5Error of x2 (τ= 0.015)

6. Conclusions

The processing of non-homogeneous terms is the key to solving the differential equations in dynamical systems and vibration engineering. In this paper, the Chebyshev’s interpolation theory is applied to the polynomial approximation of non-homogeneous terms, and an accurate numerical solution can be obtained. The method is stable, and the precision is controllable.

It can be seen from the numerical examples that when the Chebyshev interpolation nodes are used, even if a larger integration step is selected, the error of the method is far less than the error of the literature [4]. And the interpolation order selected in this paper is the same as the literature [4]. If you want to achieve higher precision requirements, you only need to select the appropriate fitting order according to the accuracy design method described in this paper. Therefore, the method can be applied to any non-homogeneous dynamical systems and vibration engineering problems.