# A homogeneous high precision direct integration based on Chebyshev interpolation

## Dongliang Chen1, Lianfu Wang2, Jindong Zhang3

1, 2, 3Harbin Engineering University, Harbin, China

2Corresponding author

Vibroengineering PROCEDIA, Vol. 22, 2019, p. 170-175. https://doi.org/10.21595/vp.2018.20444
Received 9 December 2018; accepted 16 December 2018; published 15 March 2019

Copyright © 2019 Dongliang Chen, et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
Views 70
Abstract.

Based on Chebyshev’s interpolation theory, the non-homogeneous term of the second-order linear differential equations is interpolated, and a precise integration algorithm with easy programming, high computational efficiency and precision design is realized. The method does not involve inverse operation, and does not need to additionally calculate the matrix index on the integration point, and can control the error boundary based on different precision requirements, so it has high stability and controllability. Numerical examples of periodic loads common in vibration engineering show the effectiveness of the method.

Keywords: high precision direct integration, homogeneous expansion, Chebyshev interpolation, precision design.

#### 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  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.  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.  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.  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:

(1)
$\left\{\begin{array}{l}\stackrel{˙}{v}\left(t\right)=Hv\left(t\right)+r\left(t\right),\\ v\left(0\right)={v}_{0},\end{array}\right\$

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:

(2)
$v\left(t\right)=\mathrm{exp}\left[H\bullet \left(t-{t}_{0}\right)\right]\bullet {v}_{0}.$

Let the integration step size $t-{t}_{0}=\tau$, then:

(3)
$v\left(t\right)=\mathrm{exp}\left(H\bullet \tau \right)\bullet {v}_{0}=T\bullet {v}_{0},$

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

(4)

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

(5)
$T=\mathrm{exp}\left(H\bullet \tau \right)={\left[\mathrm{e}\mathrm{x}\mathrm{p}\left(H\bullet \tau /m\right)\right]}^{m}.$

Let $m={2}^{N}$ and the original integration step size become $∆t=\tau /m$. According to the literature , if $N=$ 20, 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:

(6)
$\mathrm{exp}\left(H\bullet ∆t\right)\approx I+H\bullet ∆t+\frac{{\left(H\bullet ∆t\right)}^{2}}{2}=I+{T}_{a},$

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}$:

(7)
$\left(I+{T}_{a}\right)\left(I+{T}_{a}\right)=I+2{T}_{a}+{T}_{a}×{T}_{a}=I+{T}_{b},$

where ${T}_{b}=2{T}_{a}+{T}_{a}×{T}_{a}$. So, calculating $T$ is equivalent to executing a loop statement:

(8)

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

(9)
$T=I+{T}_{a}.$

Fig. 1. Schematic 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:

(10)

when performing interpolation point sampling, if uniform sampling is performed as described in literature , 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 $\left[a$, $b\right]$, the node coordinates are respectively:

(11)
$\left\{\begin{array}{l}{x}_{i}=\frac{b+a}{2}+\frac{b-a}{2}\mathrm{c}\mathrm{o}\mathrm{s}\frac{\left(2i-1\right)\pi }{2n},\\ i=1,2,\dots ,n.\end{array}\right\$

Then use the Newton’s divided difference formula  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:

(12)
$\frac{d}{dt}\left\{\begin{array}{c}v\left(t\right)\\ 1\\ t\\ {t}^{2}\end{array}\right\}=\left[\begin{array}{cccc}H& {r}_{0}& {r}_{1}& {r}_{2}\\ 0& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 2& 0\end{array}\right]\left\{\begin{array}{c}v\left(t\right)\\ 1\\ t\\ {t}^{2}\end{array}\right\}.$

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 ,\left({x}_{n},{y}_{n}\right)$ are fitted, and the interpolation error is:

(13)
$\left\{\begin{array}{l}f\left(x\right)-P\left(x\right)=\frac{\left(x-{x}_{1}\right)\left(x-{x}_{2}\right)\cdots \left(x-{x}_{n}\right)}{n!}{f}^{\left(n\right)}\left(c\right),\\ \mathrm{min}\left(x,{x}_{1},\dots ,{x}_{n}\right)

There is an inequality according to Chebyshev’s theory:

(14)
$\left|\left(x-{x}_{1}\right)\cdots \left(x-{x}_{n}\right)\right|\le \frac{{\left(\frac{b-a}{2}\right)}^{n}}{{2}^{n-1}},$

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

(15)
$\left|f\left(x\right)-P\left(x\right)\right|=\frac{\left|\left(x-{x}_{1}\right)\cdots \left(x-{x}_{n}\right)\right|}{n!}\left|{f}^{\left(n\right)}\left(c\right)\right|\le \frac{{\left(\frac{b-a}{2}\right)}^{n}}{n!{2}^{n-1}}\left|{f}^{\left(n\right)}\left(c\right)\right|.$

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:

(16)

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

(17)
$\left\{\begin{array}{l}{v}_{k+1}=T{v}_{k}-T\left[a\mathrm{sin}\left(\omega {t}_{k}\right)+b\mathrm{c}\mathrm{o}\mathrm{s}\left(\omega {t}_{k}\right)\right]+a\mathrm{sin}\left(\omega {t}_{k+1}\right)+b\mathrm{c}\mathrm{o}\mathrm{s}\left(\omega {t}_{k+1}\right),\\ a={\left({\omega }^{2}I+{H}^{2}\right)}^{-1}\left(\omega {r}_{2}-H{r}_{1}\right),\\ b={\left({\omega }^{2}I+{H}^{2}\right)}^{-1}\left(-\omega {r}_{1}-H{r}_{2}\right),\end{array}\right\$

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 and transforming the Eq. (16) into a first-order differential equation:

(18)
$\frac{d}{dt}\left\{\begin{array}{c}\begin{array}{c}{x}_{1}\left(t\right)\\ {x}_{2}\left(t\right)\end{array}\\ \begin{array}{c}{\stackrel{˙}{x}}_{1}\left(t\right)\\ {\stackrel{˙}{x}}_{2}\left(t\right)\end{array}\end{array}\right\}=\left[\begin{array}{cc}\begin{array}{cc}\begin{array}{c}\begin{array}{c}0\\ 0\end{array}\\ \begin{array}{c}-3\\ 2\end{array}\end{array}& \begin{array}{c}\begin{array}{c}0\\ 0\end{array}\\ \begin{array}{c}1\\ -4\end{array}\end{array}\end{array}& \begin{array}{cc}\begin{array}{c}\begin{array}{c}1\\ 0\end{array}\\ \begin{array}{c}0\\ 0\end{array}\end{array}& \begin{array}{c}\begin{array}{c}0\\ 1\end{array}\\ \begin{array}{c}0\\ 0\end{array}\end{array}\end{array}\end{array}\right]\left\{\begin{array}{c}\begin{array}{c}{x}_{1}\left(t\right)\\ {x}_{2}\left(t\right)\end{array}\\ \begin{array}{c}{\stackrel{˙}{x}}_{1}\left(t\right)\\ {\stackrel{˙}{x}}_{2}\left(t\right)\end{array}\end{array}\right\}+\left[\begin{array}{c}\begin{array}{c}0\\ 0\end{array}\\ \begin{array}{c}0\\ 10\end{array}\end{array}\right]\mathrm{sin}\omega t+\left[\begin{array}{c}\begin{array}{c}5\\ 0\end{array}\\ \begin{array}{c}0\\ 0\end{array}\end{array}\right]\mathrm{c}\mathrm{o}\mathrm{s}\omega t.$

Then use the methods in literature , literature , literature  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.

Table 1. Solution of ${x}_{1}$

 Step size Exact solution Literature  Literature  This article $\tau$ 0.049179094769134 0.048773943800022 0.049179261646032 0.049179157376206 $\text{2}\tau$ 0.093529329224097 0.090431813890316 0.093534457158153 0.093531255674386 $\text{3}\tau$ 0.128694562148165 0.119023415017619 0.128730881129231 0.128708237262728 $\text{4}\tau$ 0.151217817456480 0.130769583117008 0.151355950554894 0.151269993788613 $\text{5}\tau$ 0.158879713634259 0.124766264942132 0.159245530921237 0.159018459454587 $\text{6}\tau$ 0.150915748723294 0.103285473198170 0.151666720827195 0.151202014701130 $\text{7}\tau$ 0.128091171464306 0.071622645861488 0.129337190277955 0.128569019564925 $\text{8}\tau$ 0.092626109662066 0.037498539024419 0.094280809634630 0.093265231184200 $\text{9}\tau$ 0.047978284706641 0.010082245896349 0.049539060894824 0.048586427269383

It can be seen from the Fig. 2 and Fig. 3 that because the literature  uses a linear approximation, so a large error will occur as the integration step size increases. The error of the literature  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  and this paper are shown in Figs. 4, 5.

Table 2. Solution of ${x}_{2}$

 Step size Exact solution Literature  Literature  This article $\tau$ 0.000053759456644 0.000053148018043 0.000053802149951 0.000053775476829 $\text{2}\tau$ 0.000423725534198 0.000404520403817 0.000425056419986 0.000424225882643 $\text{3}\tau$ 0.001395020790211 0.001253656462584 0.001404688266858 0.001398666646464 $\text{4}\tau$ 0.003193642440144 0.002623526891556 0.003231859612234 0.003208116796166 $\text{5}\tau$ 0.005964208798938 0.004320682034443 0.006071292005589 0.006004981148495 $\text{6}\tau$ 0.009755671836064 0.005944066131288 0.009994373977258 0.009847119715415 $\text{7}\tau$ 0.014516394458498 0.006944525925839 0.014965191794752 0.014689486915063 $\text{8}\tau$ 0.020099074451105 0.006727765439525 0.020832484435684 0.020383841886791 $\text{9}\tau$ 0.026275033144355 0.004785022016463 0.027328318717349 0.026686325198872

Fig. 2. Error of ${x}_{1}$ Fig. 3. Error of ${x}_{2}$ Fig. 4. Error of ${x}_{1}$ ($\tau =$ 0.015) Fig. 5. Error of ${x}_{2}$ ($\tau =$ 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 . And the interpolation order selected in this paper is the same as the literature . 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.

#### Acknowledgements

This paper is funded by the International Exchange Program of Harbin Engineering University for Innovation-oriented Talents Cultivation.

1. Zhong W. X., Williams F. W. A precise time step integration method. Proceedings of the Institution of Mechanical Engineers, Part C: Journal of Mechanical Engineering Science, Vol. 208, Issue 6, 1994, p. 427-430. [Search CrossRef]
2. Lin J., Shen W., Williams F. W. Accurate high-speed computation of non-stationary random structural response. Engineering Structures, Vol. 19, Issue 7, 1997, p. 586-593. [Publisher]
3. Wang M., Au F. T. K. Assessment and improvement of precise time step integration method. Computers and Structures, Vol. 84, Issue 12, 2006, p. 779-786. [Publisher]
4. Yu X., Yuying H., Jianqiang H. A method of homogenization of high precision direct integration. Journal-Huazhong University of Science and Technology Nature Science Chinese Edition, Vol. 30, Issue 11, 2002, p. 74-76. [Search CrossRef]
5. Xiang Y., Huang Y. Y., Zeng W. G. Error analysis and accuracy design for the precise time-integration method. Chinese Journal of Computational Mechanics, Vol. 19, Issue 3, 2002, p. 276-280. [Search CrossRef]
6. Sauer T. Numerical Analysis. Second Edition, Addison-Wesley, New Jersey, 2012. [Search CrossRef]

#### Cited By

 Complexity Mingfang Chen, Kaixiang Zhang, Sen Wang, Fei Liu, Jinxin Liu, Yongxia Zhang 2020
##### JVE Journals is rebranding to Extrica  