A modified homogenized highly precise direct integration method for time-varying nonhomogeneous systems

Abstract

This article investigates a novel method for time-varying systems with an impact term; we call it a modified homogenized highly precise direct integration method. Modified homogenized highly precise direct integration method can deal with the time-varying nonhomogeneous systems effectively by direct integration with high precision. Even though it is often difficult to select the time step size of integration properly for stiff problems, modified homogenized highly precise direct integration method can effectively deal with the this problem with a large time step size. By introducing new variants twice, modified homogenized highly precise direct integration method can easily deal with the nonhomogeneous term by a novel way, inherit the advantages of highly precise direct integration method and avoid the calculation of inverse matrix. The convergency and efficiency analyses are given in this article. Several numerical simulations and an application example are presented to demonstrate the high accuracy, effectiveness and application for engineering problems of modified homogenized highly precise direct integration method.

Keywords

Time-varying system highly precise direct integration dynamical structure initial value problem

Introduction

Time-varying dynamic problems arise in many areas of structural analysis, such as vehicle bridge interactive systems, rockets or aircraft systems. The dynamic response of the structure can be analysed by the direct integration methods (the Newmark-α method,¹ Wilson-θ method² or Houbolt schemes) or transferring the system into first-order evaluation problems (the Runge–Kutta, spectral deferred correction,^3,4 time discontinuous Galerkin methods⁵). However, for some stiff problems,⁶ the time step size should be selected carefully, which makes these methods time-consuming.

Highly precise direct integration (HPD)⁷ method was proposed for homogeneous linear time-invariant dynamic systems. Due to its advantages of high precision and high efficiency, HPD had been widely applied to many problems. When applied to the linear ordinary differential equations (ODEs), HPD provided the numerical solution up to the machine precision, and it was an efficient and stable method. In order to solve the nonhomogeneous systems with high precision, expanding dimension methods^8,9 were proposed. These methods deal with the nonhomogeneous terms by different function expansions, such as Fourier, Legendre or Chebyshev expansions. However, for expanding dimension methods, there were several drawbacks including too many variants introduced, large dimension and low accuracy. In this article, we propose modified homogenized method, which can avoid these drawbacks.

However, perturbation method was presented in the study of Tan and Zhong¹⁰ recently, where the original time-varying system was decomposed into two sub-Hamiltonian systems, that is, a zero-order system and a perturbation one. A high-order multiplication perturbation method (HOMPM)¹¹ was presented, in which a low-order perturbation system was decomposed into a low-order system and a high-perturbation system by variable transformation. Each subsystem can be solved by HPD, and the solution can be determined through a series of inverse transformations. However, it is time-consuming to compute the high precise matrix in every small interval. To solve variable coefficient singularly perturbed two-point boundary value problems, HOMPM¹² was put forward. But HOMPM encounter the similar problems as HPD. Another idea to improve the accuracy is time-adaptive methods. For mass-variant systems (the flight of rockets), adaptive methods were presented in the study of Wang and colleagues,^13,14 which selected the time step size or the order of the proposed method adaptively. Methods^15,16 for periodic time-varying systems were put forward, in which the precise transferred matrix in the first period was expressed by the theory of Peano–Baker series and used in the following periods. This method was only suitable for the periodic problems. For general problems, it may cost too much computation and computer’s memory.

In this article, we propose a modified homogenized highly precise direct (MHHPD) integration method. This method is efficient for a class of time-varying nonhomogeneous systems, whose coefficients matrices are polynomials, and nonhomogeneous term is complicated. Briefly, we can summarize the method into three steps as follows:

Transfer the nonhomogeneous systems into homogeneous ones by expanding dimensions via novel methods;

Convert the transferred time-varying homogeneous systems into time-invariant systems by introducing several variants;

Utilize the standard HPD method for the transferred homogeneous time-invariant systems.

The transfer matrix is generated once and used forever, which is efficient and accurate. The accuracy of the proposed method is studied and proved by comparing with those of the others which are used commonly (e.g. RK4, Radau5,¹⁷ Newmark-α,¹ ode15s,¹⁸ Wilson-θ²).

Time-varying systems and their homogenized forms

In this section, we will provide the basic equation of structural dynamic problems with four different types of nonhomogeneous terms and their related homogenized forms. In the past decades, expanding dimension methods^8,19–22 were widely used for the nonhomogeneous dynamic systems. The nonhomogeneous loading form could be classified into three types: the polynomial loading form, the sinusoidal loading form, and the arbitrary loading form. However, the accuracy of these expanding dimension methods relates to the serials expansion caused by the truncation errors. For complicated loading terms, many dimensions need to be introduced, which is time and memory consuming. Therefore, a new type expanding dimension method is proposed in this section. For some typical loading forms, MHHPD are much more accurate.

The time-varying structural dynamic systems

The equation of motion for N-degree-of-freedom (N-DOF) time-varying dynamic systems can be expressed in a matrix form as

{\begin{matrix} M (t) X ″ + C (t) X' + K (t) X = F (t) \\ X (0) = X_{0}, X' (0) = {X'}_{0} \end{matrix}

(1)

where M(t), C(t), K(t) ∈ $R^{N \times N}$ are the time-varying mass, damping and stiffness matrices, respectively. F(t) is the loading term. $X \in R^{N}$ is a vector of physical field (e.g. displacement, temperature or velocity). Symbol prime means the derivative of time.

Let $Y = X', u = (\begin{matrix} X \\ Y \end{matrix})$ and $\tilde{F} = (\begin{matrix} 0 \\ M {(t)}^{- 1} F (t) \end{matrix})$ ; then, equation (1) can be formulated into first-order equations

u' = H (t) u + \tilde{F} (t)

(2)

where $H (t) = (\begin{matrix} 0 & I \\ - M {(t)}^{- 1} K (t) & - M {(t)}^{- 1} C (t) \end{matrix})$ , and I is an identity matrix.

The novel homogenized method

For the nonhomogeneous system (equation (2)), function expansions are the typical ways for the nonhomogeneous term. For the complicated loading forms F(t), many variants have to be introduced for high precision because of the truncation error. In this subsection, we will discuss four different kinds of complicated loading forms F(t) and their homogenized forms, which avoid introducing so many variants.

Case 1. $F (t) = e^{P_{n} (t)}$ loading form.

If p_n(t) is a polynomial vector function of degree not exceeding n. In this case, p_n(t) = p_ntⁿ + p_n–1t^n–¹+··· + p₁t + p₀, the right hand side of equation (2) has the form $F (t) = e^{P_{n} (t)}$ , then, equation (2) can be written as

{\tilde{u}}^{'} (t) : = {(\begin{matrix} u (t) \\ F (t) \end{matrix})}^{'} = \tilde{A} (t) \tilde{u} (t)

(3)

where $F (t) = e^{P_{n} (t)}, \tilde{A} (t) = (\begin{matrix} H & I \\ 0 & diag ({P'}_{n} (t)) \end{matrix}) and \tilde{u} (t) = (\begin{matrix} u (t) \\ F (t) \end{matrix})$

Case 2. $F (t) = e^{P_{n} (t)} q_{m} (t)$ loading form.

If p_n(t) is a polynomial scalar function of degree not exceeding n and q_m(t) is a polynomial vector function of degree not exceeding m, then, the nonhomogeneous system (equation (2)) can be written as

\tilde{u'} (t) : = {(\begin{matrix} u (t) \\ F (t) \end{matrix})}^{'} = \tilde{A} (t) \tilde{u} (t)

(4)

where $F (t) = e^{P_{n} (t)}, \tilde{A} (t) = (\begin{matrix} H & q_{m} (t) \\ 0 & p'_{n} (t) + 1 \end{matrix}) and \tilde{u} (t) = (\begin{matrix} u (t) \\ F (t) \end{matrix})$ .

Case 3. F(t) = p_n(t) sin(wt)+q(t)_m cos(wt) loading form.

If p_n(t) is a polynomial vector function of degree not exceeding n and q_m(t) is a polynomial vector function of degree not exceeding m, then, the nonhomogeneous system (equation (2)) can be written as

\tilde{u'} (t) : = {(\begin{matrix} u (t) \\ \sin (w t) \\ \cos (w t) \end{matrix})}^{'} = \tilde{A} (t) \tilde{u} (t)

(5)

where $\tilde{A} (t) = (\begin{matrix} H & p_{n} (t) & q_{m} (t) \\ 0 & 0 & w \\ 0 & - w & 0 \end{matrix}) and \tilde{u} (t) = (\begin{matrix} u (t) \\ \sin (w t) \\ \cos (w t) \end{matrix})$ .

Case 4. $F (t) = e^{P_{n} (t)} (q_{m} (t) \sin (w t) + r_{l} (t) \cos (w t))$ loading form.

If p_n(t) is a polynomial scalar function of degree not exceeding n, q_m(t) is a polynomial vector function of degree not exceeding m and r_l(t) is a polynomial vector function of degree not exceeding l, then, the nonhomogeneous system (equation (2)) can be written as

\tilde{u'} (t) : = {(\begin{matrix} u (t) \\ P_{1} (t) \\ P_{0} (t) \end{matrix})}^{'} = \tilde{A} (t) \tilde{u} (t)

(6)

where $P_{1} (t) = e^{P_{n} (t)} \sin (w t), P_{0} (t) = e^{P_{n} (t)} \cos (w t), \tilde{A} (t) = (\begin{matrix} H & q_{m} (t) & r_{l} (t) \\ 0 & P'_{n} (t) & w \\ 0 & - w & P'_{n} (t) \end{matrix}) and \tilde{u} (t) = (\begin{matrix} u (t) \\ P_{1} (t) \\ P_{0} (t) \end{matrix})$

From equations (3)–(6), we can summarize that the nonhomogeneous systems can be homogenized into the same form, as follows

{\begin{matrix} \tilde{u'} (t) = \tilde{A} (t) \tilde{u} (t) \\ \tilde{u} |_{t = t_{0}} = \tilde{u} (t_{0}) \end{matrix}

(7)

where $\tilde{u}$ and $\tilde{A}$ are the expanding forms of u and H, respectively. The extra dimensions depend on different forms of F(t). The coefficients matrix $\tilde{A}$ is a polynomial matrix, and the form of $\tilde{A}$ depends on the time-varying terms of equation (1). MHHPD is suitable for the problems whose coefficients based on polynomials

\tilde{A} = \sum_{i = 0}^{r} {\tilde{A}}_{i} t^{i}

MHHPD

In this section, we introduce MHHPD more in detail and analyse the convergence of this method. The nonhomogeneous system (equation (2)) with coefficient A can be transferred into the homogenized one (equation (7)) whose coefficient is time-varying $\tilde{A}$ . For the convenience of convergence analysis, we have to transfer the homogenized system into an unit time-varying one.

Definition 1

For 0 < ρ < 1 and $A (s) = \sum_{i = 0}^{r} A_{i} s^{i}$ , if

{\begin{matrix} u' (s) = \frac{1}{ρ} A (s) u (s) \\ u |_{s = 0} = u (0) \end{matrix} s ϵ [0, ρ]

(8)

Then, equation (8) is called unit time-varying system, A_i is an n × n time-invariant coefficient matrix, and A(s) is time-varying coefficient matrix based on s, 0 < s < ρ < 1. The system (equation (7)) can be transformed into the unit time-varying system by the following mapping

s = ρ \frac{t - t_{0}}{t_{f} - t_{0}}

(9)

Then, the time interval [t₀, t_f] can be mapped into [0, ρ], with the initial value s(t₀) = 0. According to equation (9), we can get the related u(s) and A(s) as well

\begin{matrix} u (s) = \tilde{u} (t) = \tilde{u} (t_{0} + \frac{1}{ρ} (t_{f} - t_{0}) s) \\ A (s) = \tilde{A} (t) = \tilde{A} (t_{0} + \frac{1}{ρ} (t_{f} - t_{0}) s) \end{matrix}

We can make the conclusion that system equations (7) and (8) are equivalent.

Expanding dimensions for the time-varying systems

This subsection introduces how to deal with the time-varying coefficient $\tilde{A} (t)$ . Some new variables should be introduced as well.

Let

\begin{matrix} {\tilde{u}}_{0} {(s)}_{n \times 1} = \tilde{u} {(s)}_{n \times 1} \\ {\tilde{u}}_{i} (s) = s^{i} {\tilde{u}}_{0} (s), (i = 1, 2, \dots, m + r; m > r) \end{matrix}

(10)

and

\hat{U} {(s)}_{n (m + r + 1) \times 1} = {({\tilde{u}}_{0}^{T} {\tilde{u}}_{1}^{T} \dots {\tilde{u}}_{r}^{T} \dots {\tilde{u}}_{m}^{T} {\tilde{u}}_{m + 1}^{T} \dots {\tilde{u}}_{m + r}^{T})}^{T}

(11)

where, r is the degree of $\tilde{A} (t)$ and m+r > M (M is a non-negative integer). Moreover

\begin{matrix} \tilde{u}'_{i} (s) = i s^{i - 1} u_{0} (s) + s^{i} u'_{0} (s) = i s^{i - 1} u_{0} (s) + A (s) (s^{i} u_{0} (s)) \\ = i I_{n} u_{i - 1} (s) + A (s) u_{i} (s), (i = 1, 2, \dots, m) \end{matrix}

(12)

\begin{matrix} {\tilde{u'}}_{i} (s) = i s^{i - 1} u_{0} (s) + s^{i} u'_{0} (s) \\ = i I_{n} u_{i - 1} (s) + s^{i} I_{n} u'_{0} (s), (i = m + 1, \dots, m + r) \end{matrix}

(13)

Hence, a combination of equations (8)–(13) leads to

{\begin{matrix} W \hat{U'} (s) = \hat{A} \hat{U} (s) \\ \hat{U} |_{s = 0} = {(u {(0)}^{T}, 0, \dots, 0)}^{T} \end{matrix} s ϵ [0, ρ]

(14)

where

Then, the time-varying coefficient matrix A(s) is transferred into time-invariant one $\hat{A}$ .

Convergence

In this subsection, we will study the convergence of MHHPD. To this end, we give the following.

Definition 2

If $\overset{\lor}{A} = \hat{A}$ is an time-invariant matrix, then

{\begin{matrix} {\overset{\lor}{U}}^{'} (s) = \overset{\lor}{A} \overset{\lor}{U} (s), \\ \overset{\lor}{U} |_{s = 0} = {(u {(0)}^{T}, 0, \dots, 0)}^{T}, \end{matrix} s ϵ [0, ρ],

(15)

is called time-invariant homogeneous system, where $\overset{\lor}{U} = \hat{U}$ , $\overset{\lor}{U} |_{s = 0} = \hat{U} |_{s = 0}$ .

In the views of some studies,^9,21 equations (2) and (7) are equivalent for the expanding dimension method.

Theorem 1

When ε(t) → 0, the solution of

\tilde{u'} = (A (t) + ε (t)) \tilde{u} + F (t)

(16)

converges to the solution

u' = (A (t)) u + F (t)

(17)

Proof

By equations (16) and (17), we obtain

{\begin{matrix} {(Δ u)}^{'} = A (t) Δ u + ε (t) \tilde{u} \\ (Δ u) (0) = 0 \end{matrix}

(18)

where $Δ u = \tilde{u} - u$ .

Let $ε (t) \tilde{u} = g (t)$ , then

{\begin{matrix} {(Δ u)}^{'} = A (t) Δ u + g (t) \\ (Δ u) (0) = 0 \end{matrix}

(19)

Hence

‖ Δ u ‖_{\infty} \leq C ‖ g (t) ‖_{\infty} \leq \tilde{C} ‖ ε (t) ‖_{\infty}

So, as ε(t) → 0, the solution of equation (16) converges to the solution of equation (17).

Theorem 2

When m → +∞, the solution of equation (14) converges to the solution of equation (15).

Proof

A combination of equations (14) and (15) leads to

{\begin{matrix} \hat{U'} (s) - {\overset{\lor}{U}}^{'} (s) = W^{- 1} \hat{A} \hat{U} (s) - \overset{\lor}{A} \overset{\lor}{U} (s) \\ (\hat{U} - \overset{\lor}{U}) |_{s = 0} = {(0, 0, \dots, 0)}^{T} \end{matrix} s ϵ [0, ρ]

Furthermore, we partition matrices W and $\hat{A}$ , let

\begin{array}{l} W = (\begin{matrix} I_{n (m + 1) \times n (m + 1)} & 0_{n (m + 1) \times n r} \\ - W_{ε} & I_{n r \times n r} \end{matrix}), \\ \hat{A} = (\begin{matrix} B_{n (m + 1) \times n (m + 1)} & C_{n (m + 1) \times n r} \\ 0_{n r \times n (m + 1)} & D_{n r \times n r} \end{matrix}) \end{array}

where

W_{ε} = {(\begin{matrix} \begin{matrix} \begin{matrix} - s^{m + 1} I & 0 \end{matrix} & \begin{matrix} \dots & 0 \end{matrix} \end{matrix} \\ \begin{matrix} \begin{matrix} \begin{matrix} ⋮ & ⋮ \end{matrix} \\ \begin{matrix} - s^{m + r} I & 0 \end{matrix} \end{matrix} & \begin{matrix} \begin{matrix} ⋮ & ⋮ \end{matrix} \\ \begin{matrix} \dots & 0 \end{matrix} \end{matrix} \end{matrix} \end{matrix})}_{n r \times n (m + 1)}

It follows that

W^{- 1} = (\begin{matrix} I_{n (m + 1) \times n (m + 1)} & 0_{n (m + 1) \times n r} \\ W_{ε} & I_{n r \times n r} \end{matrix})

Moreover

W^{- 1} \hat{A} = (\begin{matrix} B & C \\ W_{ε} B & W_{ε} C + D \end{matrix})

Therefore

\begin{matrix} \hat{U'} (s) - {\overset{\lor}{U}}^{'} (s) = W^{- 1} \hat{A} \hat{U} (s) - \hat{A} \overset{\lor}{U} (s) \\ = (\begin{matrix} B & C \\ W_{ε} B & W_{ε} C + D \end{matrix}) \hat{U} (s) - (\begin{matrix} B & C \\ 0 & D \end{matrix}) \overset{\lor}{U} (s) \\ = (\begin{matrix} B & C \\ 0 & D \end{matrix}) (\hat{U} (s) - \overset{\lor}{U} (s)) + (\begin{matrix} 0 & 0 \\ W_{ε} B & W_{ε} C \end{matrix}) \hat{U} (s) \end{matrix}

This, together with equation (8), leads to

‖ (\begin{matrix} 0 & 0 \\ W_{ε} B & W_{ε} C \end{matrix}) ‖_{1} = lim_{m \to + \infty} (s^{m + 1} {‖ {\hat{A}}_{i} ‖}_{1}) \to 0

Hence, we obtain

{\hat{U}}^{'} (s) - {\overset{\lor}{U}}^{'} (s) = (\begin{matrix} B & C \\ 0 & D \end{matrix}) (\hat{U} (s) - \overset{\lor}{U} (s))

with the initial conditions

(\hat{U} (s) - \overset{\lor}{U} (s)) |_{s = 0} = {(0, 0, \dots, 0)}^{T}

Finally, all the above yields

\hat{U} (s) \to \overset{\lor}{U} (s)

HPD method for the transferred systems

Equation (15) is a homogenized time-invariant system, standard HPD can be used. The solution of equation (15) can be written as

u (t) = \exp (\hat{A} \times t) u (0)

If τ is taken as step size of the time, that is, $τ = t_{k + 1} - t_{k} (k = 0, 1, 2, \dots)$ , then

u (t_{k + 1}) = T (τ) u (t_{k})

where

T (τ) = \exp (\hat{A} \times τ)

As $u (t_{k})$ is known from the calculation of the previous time step, the computation of matrix T(τ) is crucial for determining $u (t_{k + 1})$ .

In order to compute T(τ) accurately, it is desirable to transform

T (τ) = \exp (\hat{A} \times τ) = {(\exp (\hat{A} \times d τ))}^{2^{N}}

(20)

where $d τ = τ / 2^{N}$ .

Furthermore, the precision parameter N should be selected previously (e.g. N = 20). Then

T (τ) = {(I + d τ \times \hat{A})}^{2^{N}} = : {(I + T_{α 0})}^{2^{N}}

Furthermore

\begin{matrix} \begin{matrix} {(I + T_{α 0})}^{2} = I + 2 T_{α 0} + T_{α 0} \times T_{α 0} = : I + T_{α 1} \\ {(I + T_{α 1})}^{2} = I + 2 T_{α 1} + T_{α 1} \times T_{α 1} = : I + T_{α 2} \end{matrix} \\ \begin{matrix} \dots \dots \dots \\ {(I + T_{α, N - 1})}^{2} = I + 2 T_{α, N - 1} + T_{α, N - 1} \times T_{α, N - 1} = : I + T_{α, N} \end{matrix} \end{matrix}

It is clear that

I + T_{α, N} = {(I + T_{α, N - 1})}^{2} = \dots = {(I + T_{α, 0})}^{2^{N}} = T (τ)

The procedure for computing T(τ) is described in Algorithm 1 HPD.

Algorithm 1.

HPD.

d τ = τ / 2^{N}

T_{α} = d τ \times \hat{A} + \frac{1}{2} {(d τ \times \hat{A})}^{2}

R e p e a t

N = 0
N = N + 1

T_{α} = T_{α} \times T_{α} + 2 T_{α}

Unti lN = 30, % for example

T = I + T_{α}

Finally, we give the procedure of MHHPD in Algorithm 2 MHHPD.

Algorithm 2.

MHHPD.

Step 1: Homogenized the loading form
Switch

F (t)

case 1

simple form^[21]

case 2 e^{p_{n} (t)}

: refer to equation (3)

case 3 e^{p_{n} (t)} q_{m} (t)

: refer to equation (4)

case 4 p_{n} (t) \sin (w t) + q_{m} (t) \cos (w t)

: refer to equation (5)

case 5 e^{p_{n} (t)} (q_{m} (t) \sin (w t) + r_{l} (t) \cos (w t))

: refer to equation (6)
End Switch
Step 2: Transfer the time-varying systems into time-invariant ones

m = 20, r = 2

, % for example
Introduce some variables, refer to equation (10)
Get the time-invariant matrix

\hat{A}

, refer to equation (14)
Step 3: Algorithm 1 HPD

Efficiency analysis

For the dynamical problems, especially, the stiff problems, we have to select a small time step size for stability. There is no doubt, small time step size costs lots of computation. MHHPD method inherits the advantages of HPD method and can calculate the stiff problem by large time step. Even though MHHPD is a bit complicated before using the standard HPD method, it is worthy for us compared to the time cost with different methods. Therefore, in this subsection, we roughly estimate the cost time of MHHPD and compare the cost with traditional one-step method (e.g. Newmark-β method). To begin with, we assume that the cost of each step is identical for both MHHPD and traditional one.

First, we assume that the high precise parameter is N, and the time step size is dt. We take the $OneSte p_{t}$ as the cost time of each time step. For Algorithm 2, we have to calculate N times with sub-time step size $τ = d t / 2^{N}$ , and we know the time cost for every time interval is

M H H P D_{t} = N * OneSte p_{t}

However, for traditional methods, 2^N sub-steps have to be calculated for equivalent accuracy within every time interval

T r a d i t i o n_{t} = 2^{N} * OneSte p_{t}

Based on the analysis hereinbefore, the time cost of the traditional method is about $2^{N} / N$ times as much as MHHPD.

Solving process

Let t = t_k₊₁. The numerical solution of equation (15) is

u (t_{k + 1}) = T (τ) u (t_{k}) = \exp (\hat{A} \times τ) u (t_{k})

Using Algorithm 1 HPD, the solution u(t_k₊₁) can be obtained with high precision. Here, the step size τ can cover an arbitrary segment. By selecting a proper parameter N, though we select a large step size, the precise solution will be obtained. Since equation (15) is a linear time-invariant system, the transfer matrix is computed once and can be used forever, which saves a lot computation.

Finally, the unit interval needs to be transformed into the original interval by an inverse transform

\tilde{u} (t) = u (s) = u (ρ \frac{t - t_{0}}{t_{f} - t_{0}})

Numerical simulations

In this section, we present several numerical results to illustrate the advantages of MHHPD.

A first-order time-varying nonhomogeneous system

We show the steps and efficiency of MHHPD by a simple example for a better understanding

{\begin{matrix} u' (t) + t u (t) = \frac{1}{3} t^{4} + \frac{1}{2} t^{3} + 2 t^{2} + 2 t + 1 \\ u (0) = 1 \end{matrix}

(21)

The exact solution is given by

u_{e} = \frac{1}{3} t^{3} + \frac{1}{2} t^{2} + t + 1

We solve equation (21) using MHHPD method and fourth-order Runge–Kutta (RK4²³) method with different time step sizes.

We present several necessary steps of MHHPD in detail as follows:

Homogenize the systems. By expanding dimension method, equation (21) can be homogenized as follows

\tilde{u'} (t) = \tilde{A} (t) \tilde{u} (t) = (a_{0} + a_{1} t) \tilde{u} (t)

where

\begin{matrix} \tilde{u} (t) = (\begin{matrix} \begin{matrix} u (t) \\ t^{4} \\ t^{3} \end{matrix} \\ \begin{matrix} t^{2} \\ t \\ 1 \end{matrix} \end{matrix}), a_{1} = (\begin{matrix} \begin{matrix} 0 & \frac{1}{3} & \frac{1}{2} \\ 0 & 0 & 4 \\ 0 & 0 & 0 \end{matrix} & \begin{matrix} 2 & 2 & 1 \\ 0 & 0 & 0 \\ 3 & 0 & 0 \end{matrix} \\ \begin{matrix} 0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{matrix} & \begin{matrix} 0 & 2 & 0 \\ 0 & 0 & 1 \\ 0 & 0 & 0 \end{matrix} \end{matrix}), \\ a_{1} = (\begin{matrix} \begin{matrix} - 1 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{matrix} & \begin{matrix} 0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{matrix} \\ \begin{matrix} 0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{matrix} & \begin{matrix} 0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{matrix} \end{matrix}) \end{matrix}

Transfer the time-varying system into time-invariant one. As the highest degree of the polynomial coefficient $\tilde{A} (t)$ is r = 1, then the transferred time-invariant system can be expressed as

\overset{\lor}{u}' (s) = \overset{\lor}{A} \overset{\lor}{u} (s)

where

Utilize HPD method. Algorithm 1 can be used for the high precision.

In Table 1, the comparison of the numerical errors between RK4 and MHHPD under the same time step size is given. It can be observed that MHHPD is much more accurate than RK4. It is also shown in Table 1 that time step size h has little influence on the numerical error of MHHPD. If h is fixed, as the expand dimension parameter m increases, the error decreases. However, when m > 10, m has little influence on the accuracy.

Table 1.

Simulation 5.1 – comparison of errors with MHHPD (N = 20) and RK4.

Scheme	Parameter	Step size	0.2	0.4	0.6	0.8
rk4	−	0.2	2.6293e − 06	2.6293e − 06	2.6293e − 06	2.6293e − 06
MHHPD	$m = 4$	0.2	6.8168e − 14	1.5443e − 10	1.4744e − 08	3.8344e − 07
	$m = 6$		0.0000e+00	3.9968e − 15	2.4396e − 12	1.9998e − 10
	$m = 8$		0.0000e+00	−2.2204e − 16	−2.2204e − 16	6.3505e − 14
	$m = 10$		0.0000e+00	−2.2204e − 16	−4.4409e − 16	0.0000e+00
	$m = 15$		0.0000e+00	−2.2204e − 16	−4.4409e − 16	0.0000e+00
RK4		0.1	1.6546e − 07	4.1371e − 07	8.0179e − 07	1.4102e − 06
MHHPD	$m = 4$	0.1	6.8390e − 14	1.5443e − 10	1.4744e − 08	3.8344e − 07
	$m = 6$		4.4409e − 16	5.3291e − 15	2.4403e − 12	1.9998e − 10
	$m = 8$		4.4409e − 16	4.4409e − 16	0.0000e+00	6.0840e − 14
	$m = 10$		4.4409e − 16	4.4409e − 16	−4.4409e − 16	−1.7764e − 15
	$m = 15$		4.4409e − 16	4.4409e − 16	−4.4409e − 16	−1.7764e − 15
RK4	−	0.05	1.0352e − 08	2.5811e − 08	4.9805e − 08	8.7137e − 08
MHHPD	$m = 4$	0.05	6.7724e − 14	1.5443e − 10	1.4744e − 08	3.8344e − 07
	$m = 6$		2.2204e − 16	4.6629e − 15	2.4396e − 12	1.9999e − 10
	$m = 8$		2.2204e − 16	−2.2204e − 16	−2.2204e − 16	6.5725e − 14
	$m = 10$		2.2204e − 16	−2.2204e − 16	−6.6613e − 16	2.2204e − 15
	$m = 15$		2.2204e − 16	−2.2204e − 16	−6.6613e − 16	2.2204e − 15

MHHPD: modified homogenized highly precise direct integration method.

In Table 2, the influence of dimensional expanding parameter m and precision parameter N is discussed. For time-invariant systems, the selection of precision parameter N was discussed in Tan et al.²⁴ For the time-varying system (equation (21)), N has little influence on the accuracy. If m < 10, as m increases the error decreases; but when m > 10, the errors no longer decrease, as m increases, so there is no need for a bigger m.

Table 2.

Simulation 5.1 – comparison of errors with MHHPD (h = 0.1) with different dimensional expanding m and different parameter N.

Scheme	N	m	0.2	0.4	0.6	0.8	1
mhhpd	1	51015	8.8818e − 16	−4.5253e − 13	−1.5576e − 10	−8.2599e − 09	−1.7323e − 07
			8.8818e − 16	4.4409e − 16	6.6613e − 16	8.8818e − 16	1.7764e − 15
			8.8818e − 16	4.4409e − 16	6.6613e − 16	8.8818e − 16	4.4409e − 16
MHHPD	2	51015	−6.6613e − 16	−8.3999e − 13	−1.9517e − 10	−9.2179e − 09	−1.8429e − 07
			−6.6613e − 16	−1.3323e − 15	−1.7764e − 15	−4.4409e − 16	1.3323e − 15
			−6.6613e − 16	−1.3323e − 15	−1.7764e − 15	−4.4409e − 16	0.0000e+00
MHHPD	4	51015	0.0000e+00	−9.4569e − 13	−2.0353e − 10	−9.3989e − 09	−1.8625e − 07
			0.0000e+00	1.1102e − 15	2.2204e − 16	4.4409e − 16	3.5527e − 15
			0.0000e+00	1.1102e − 15	2.2204e − 16	4.4409e − 16	1.3323e − 15
MHHPD	6	51015	2.2204e − 16	−9.4835e − 13	−2.0369e − 10	−9.4023e − 09	−1.8628e − 07
			4.4409e − 16	1.5543e − 15	1.5543e − 15	0.0000e+00	3.1086e − 15
			4.4409e − 16	1.5543e − 15	1.5543e − 15	0.0000e+00	8.8818e − 16
MHHPD	8	51015	−4.4409e − 16	−9.4968e − 13	−2.0370e − 10	−9.4023e − 09	−1.8628e − 07
			−4.4409e − 16	0.0000e+00	−4.4409e − 16	0.0000e+00	4.8850e − 15
			−4.4409e − 16	0.0000e+00	−4.4409e − 16	0.0000e+00	2.2204e − 15

MHHPD: modified homogenized highly precise direct integration method.

A time-varying nonhomogeneous system with multi-frequency

For the multi-frequency systems, the canonical equations are stiff. So, the time step size has to be selected carefully.⁶ Since MHHPD is insensitive to the time step size, it is convenient for this problem.

Consider the system

{\begin{matrix} u'_{1} + (1 + a + t) u_{2} - a u_{3} + b u_{4} - b u_{5} = F 1 \\ \begin{matrix} u'_{2} + (t + 1) u_{2} - a u_{3} = F 2 \\ u'_{3} + a u_{2} + (t) u_{3} = F 3 \\ u'_{4} + (t + 1) u_{4} - b u_{5} = F 4 \\ u'_{5} + b u_{4} + (t) u_{5} = F 5 \end{matrix} \end{matrix}

(22)

where

{\begin{matrix} \begin{matrix} F 1 = \cos (at) + \cos (bt) + \sin (at) + \sin (bt) + t \sin (at) (t + 1) \\ F 2 = \sin (at) + t \sin (at) (t + 1) \end{matrix} \\ \begin{matrix} F 3 = \cos (at) + t^{2} \cos (at) \\ F 4 = \sin (bt) + t \sin (bt) (t + 1) \\ F 5 = \cos (bt) + t^{2} \cos (bt) \end{matrix} \end{matrix}

The exact solution is given by

\begin{matrix} u_{1} = t \cos (at) + t \cos (bt) + t \sin (at) + t \sin (bt) \\ u_{2} = t \sin (at), u_{3} = t \cos (at) \\ u_{4} = t \sin (bt), u_{4} = t \cos (bt) \end{matrix}

Equation (22) is a multi-frequency system, and the parameters a and b represent the different frequencies. We take a = 10 and b = 10,000 as an example. Few numerical methods can solve this system with big time step size. For a better comparison, we apply MHHPD (N = 35; m = 40) with time step size h = 0.01, time-adaptive methods (Radau5¹⁷ and ode15s¹⁸) and RK4 with time step size h = 10⁻⁶ for this problem.

For the multi-frequency problem, traditional methods can simulate the low frequency very well; however, for the high frequency, things changed, the time step size has to be small enough for a better solution. In our simulation, the solutions of MHHPD and traditional methods are equivalent in the aspect of low-frequency part. In Figures 1 –3, the numerical solutions by mentioned methods are presented for the high-frequency part. The right one is the exact solutions, and the left ones are the numerical solution. Radau5 is an adaptive method, and we choose the relative error and absolute error as 10⁻⁶, and the number of computational points is 6132. Similarly, ode15s is an adaptive method; the relative error and the absolute error are the same as Radau5, and the number of computational points is 119,707. RK4 is an explicit method, which needs tiny time step size h = 10⁻⁶ for a good solution. The time step size of MHHPD is relatively big h = 0.01. As shown in Figures 1 –3, although the time step size of MHHPD is big, MHHPD is more accurate than traditional methods.

Figure 1.

Simulation 5.2 – the numerical solution of u₁ obtained by compared methods: (a) the exact solution and (b) numerical solutions.

Figure 2.

Simulation 5.2 – the numerical solution of u₄ obtained by compared methods: (a) the exact solution and (b) numerical solutions.

Figure 3.

Simulation 5.2 – the numerical solution of u₅ obtained by compared methods: (a) the exact solution and (b) numerical solutions.

In Table 3, the errors and CPU time are presented. MHHPD is an explicit method and has fewer computational costs and is much more efficient than the traditional methods.

Table 3.

Simulation 5.2 – comparison of the errors by different methods.

Scheme	MHHPD	Radau5	ode15s	RK4
$\underset{0 \leq t \leq T}{m a x} ∥ u_{e} - u_{h} ∥_{\infty}$	3.4516e − 05	1.9755e − 04	3.8303e − 02	1.6272e+00
CPU time	1.6114e − 03	7.2053e+00	1.2332e+01	1.6412e+01
Total steps	160	6132	119,707	1,600,000

MHHPD: modified homogenized highly precise direct integration method.

A nonhomogeneous Hermite equation

Consider the nonhomogeneous Hermite equation

{\begin{matrix} u ″ - (1 + 2 t) u' + (1 + t) u = F \\ u |_{t = 0} = 0, u' |_{t = 0} = e \end{matrix}

(23)

We consider a complicate form of F

\begin{matrix} F = t^{3} (2 e^{1 - t} (50 \cos (50 t) - 50 \sin (50 t)) - 3 e^{1 - t} (\cos (50 t) + \sin (50 t))) \\ + t^{2} e^{1 - t} (50 \cos (50 t) - 50 \sin (50 t)) + t^{2} e^{1 - t} (2500 \cos (50 t) + 2500 \sin (50 t)) \\ + 4 t^{2} e^{1 - t} (\cos (50 t) + \sin (50 t)) - 7 t e^{1 - t} (50 \cos (50 t) - 50 \sin (50 t)) \\ - t e^{1 - t} (2500 \cos (50 t) + 2500 \sin (50 t)) + 7 t e^{1 - t} (\cos (50 t) + \sin (50 t)) \\ + 2 e^{1 - t} (50 \cos (50 t) - 50 \sin (50 t) - 5 e^{1 - t} (\cos (50 t) + \sin (50 t))) \end{matrix}

Equation (23) has the exact solution

u (t) = e^{1 - t} (t - t^{2}) (\cos (50 t) + \sin (50 t))

Taylor or Chebyshev series expansions of F(t) are recommended in the literature.^20–22 Because of truncation errors, when using the recommended methods, we have to introduce many dimensions, which unavoidably cost memories. MHHPD can avoid the large dimension and the truncation error.

This simulation has been carried out with different time integration schemes: MHHPD (N = 15;m = 30), fourth-order Runge–Kutta, Newmark method (with α = 0.25; β = 0.5),¹ Wilson method (with θ = 1.4),² Radau5 and ode15s. RK4 method, Newmark method (with α = 0.25; β = 0.5) and Wilson method (with θ = 1.4) apply the same time step size h = 10⁻⁶. Radau5 and ode15s are time-adaptive methods; we choose the relative error and absolute error as 10⁻⁸.

MHHPD is an explicit method; the transfer matrix T is calculated once and used forever, so it is efficient on this problems. In Table 4, the numerical errors and CPU time (in s) of different schemes are compared. It can be observed that the cost of CPU time of MHHPD method is much less than other methods. In addition to this, MHHPD method provides much more accurate numerical results.

Table 4.

Simulation 5.3 – the errors and CPU time by different methods.

Scheme	MHHPD	Radau5	ode15s	Newmark	Wilson	RK4
$\underset{0 \leq t \leq T}{m a x} ∥ u_{e} - u_{h} ∥_{\infty}$	3.198e − 10	1.310e − 08	1.482e − 06	1.650e − 02	1.642e − 02	2.012e − 02
$\underset{0 \leq t \leq T}{m a x} ∥ u'_{e} - u'_{h} ∥_{\infty}$	1.316e − 08	3.063e − 08	4.439e − 06	6.114e − 01	6.110e − 01	9.723e − 01
CPU time	1.320e − 03	3.091e − 01	3.640e − 01	1.122e − 01	3.432e − 01	2.745e+00

MHHPD: modified homogenized highly precise direct integration method.

The precision parameter N and dimensional expanding parameter m play important roles for the computation. The selection of N for the HPD method was discussed.²⁴ The influence of m and N is presented in Figure 4. As N increases, the sub-step τ = h/2^N and the errors decrease. However, as shown in Figure 4, when N > 15, the increase in N has little influence to the accuracy. Meanwhile, the parameter m, namely, the expanding dimensions, relates to the accuracy and has the following regularities: when m < 10, the bigger the m, the higher the accuracy, and while m ∈ (10, 20), the errors are the lowest; however, when m > 20, as m increases, the accuracy decreases, which is mainly caused by the computer’s round-off error.²⁵ Normally, the selection N, m ∈ [10, 20] is recommended for most problems. Actually, the bigger the m, the more the computer’s memory cost. For large systems, for example, the spatial discrete partial differential equation (PDE) with thousands of DOFs, the cost of computer’s memory increases tremendously.

Figure 4.

Simulation 5.3 – the errors by different parameters N and m with time step h = 0.01: (a) the displacement and (b) the velocity.

Engineering application

Dynamic system is one of the most important systems; if the coefficient matrix is fixed, the time-varying systems are decayed into time-invariant ones. For convenience, we take the time-invariant one as an example to show its application in engineering. We consider the Pratt truss, which was invented by Caleb and Thomas Pratt in 1844. The cross-sectional area, nodes, elements and length are shown in Figure 5, and the Young’s modulus E = 1e9 Pa, density ρ = 2700 kg/m³. The constant y-direction point force f ≡ –10 N is acted on point 9 with zero initial displacement and derivative. For getting the mass and stiffness matrices, we use finite element method (space bar element) to discretize this truss model. Moreover, we let the Rayleigh damping coefficients α = 2e − 2, β = 3e − 3.

Figure 5.

Example 5.4 – the Pratt truss model.

We set T = 2 s, N = 25, in the computation (cf. Algorithm 2). Owing to no exact solution for this complex structure problem, we compute a reference solution with a sufficiently small uniform time step size. Then, we compute the numerical solution using Algorithm 2, and we compare to Wilson θ (θ = 1.4) and Newmark β method (α = 0.25, β = 0.5) with uniform time step size k = 1/N^∗, which is shown in Figure 6. In Figure 7, we show the scaling deformation and the deformation value at time T = 2 s. From Figures 6 and 7, we can see that Algorithm 2 is very efficient in the Pratt truss problem. Table 5 shows that MHHPD costs less than traditional methods.

Figure 6.

Example 5.4 – numerical solutions with different method at Point 9.

Figure 7.

Example 5.4 – the scaling deformation of the Pratt truss at T = 2 s.

Table 5.

Simulation 5.3 – the time cost by different methods.

Scheme	MHHPD	ode15s	RK4
Time cost	0.3126	0.3983	0.0016
Steps	1e4	1e4	1e2
Time cost of each step	3.126e − 5	3.083e − 5	1.6e − 5

Conclusion

In this article, a MHHPD is designed for solving time-varying nonhomogeneous systems. The time-varying problems with nonhomogeneous terms cannot be implemented with high precision; however, by introducing new variables twice, we can easily deal with this problem. MHHPD inherits the advantages of HPD method, which is an explicit method. For some problems, even though some extra variables are introduced and the dimensions are expanded, the proposed method is much more effective than some traditional methods.

In subsection ‘Convergence’, the convergency is given and subsection ‘Efficiency analysis’ shows the efficiency of MHHPD method. The influence of parameters is discussed in ‘Numerical simulations’, ‘A first-order time-varying nonhomogeneous system’ and ‘A nonhomogeneous Hermite equation’; the effectiveness for stiff problem ‘A time-varying nonhomogeneous system with multi-frequency’ and ‘A nonhomogeneous Hermite equation’ is demonstrated, and the application in engineering is also shown in ‘Engineering application’.

Numerical results show that the proposed method indeed provide solutions with high precision.

Footnotes

Handling Editor: Chenguang Yang

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work was supported by the NSFC (grant no. 50876066), SSPU fund project (no.EGD15XQD14), Shanghai University Youth Teacher Training Program (no. ZZZZEGD15007) and SSPU Applied Mathematics key discipline (no. XXKPY1604).

References

Newmark

. A method of computation for structural dynamics. J Eng Mech Div 1959; 85: 67–94.

Wilson

. A computer program for the dynamic stress analysis of underground structures. Berkeley, CA: Structural Engineering Lab, University of California, 1968.

Dutt

Greengard

Rokhlin

. Spectral deferred correction methods for ordinary differential equations. BIT 2000; 40: 241–266.

Huang

Jia

Minion

. Accelerating the convergence of spectral deferred correction methods. J Comput Phys 2006; 214: 633–656.

Cockburn

Karniadakis

Shu

. The development of discontinuous Galerkin methods. In: Cockburn

Karniadakis

Shu

(eds) Discontinuous Galerkin methods. Berlin: Springer, 2000, pp. 3–50.

Cohen

Hairer

Lubich

. Numerical energy conservation for multi-frequency oscillatory differential equations. BIT 2005; 45: 287–305.

Zhong

. On precise time-integration method for structural dynamics. J Dalian Univ Technol 1994; 34: 131–136.

Chen

Zhang

et al . Precise time-integration method with dimensional expanding for structural dynamic equations. AIAA J 2001; 39: 2394–2418.

Zhou

Wang

Jia

. A homogenized high precise direct integration based on Taylor serials. J Shanghai Jiaotong Univ 2001; 35: 1916–1919.

10.

Tan

Zhong

. Numerical solution of differential Riccati equation with variable coefficients via symplectic conservative perturbation method. J Dalian Univ Technol 2006; 46: S7–S13.

11.

Lan

et al . The high order multiplication perturbation method for time-varying dynamic system. Sci Sin Phys 2012; 42: 185.

12.

Zhang

Huang

. Coupling of high order multiplication perturbation method and reduction method for variable coefficient singular perturbation problems. Appl Math Mech 2014; 35: 97–104.

13.

Wang

Xing

. Self-adaptive Newmark method of variable-mass beam dynamic system. J Beijing Univ Aeronaut Astronaut 2014; 40: 829–833.

14.

Xing

Wang

Pan

et al . Self-adaptive Lie series method for mass-variable systems. Sci Sin Technol 2014; 44: 532–536.

15.

Zhou

Jiang

. A high precision of integration method for periodically time-varying linear system. J Shanghai Jiaotong Univ 2005; 39: 1016–1019.

16.

Peng

. Solving time-varying periodic coefficient Riccati differential equations via Fourier series and precise integration method. Chinese J Comput Mech 2009; 26: 772–777.

17.

Wanner

Hairer

. Solving ordinary differential equations II. Stiff and differential-algebraic problems. 2nd ed. Berlin: Springer 1991.

18.

Shampine

Reichelt

. The MATLAB ode suite. SIAM J Sci Comput 1997; 18: 1–22.

19.

Lin

Shen

Williams

. A high precision direct integration scheme for structures subjected to transient dynamic loading. Comput Struct 1995; 56: 113–120.

20.

Zhou

Luo

et al . A homogenized high precision direct integration based on Chebyshev polynomial series. J Donghua Univ 2006; 32: 46–49.

21.

Wang

Tian

Zhou

. Homogenized high precision direct integration scheme and its applications in engineering. Int J Numer Meth Bio 2002; 18: 429–439.

22.

Shi

Zhou

. A homogeneous high precise direct integration based on legendre polynomial series. Chinese J Comput Mech 2005; 22: 335–338.

23.

Burden

Faires

. Numerical analysis. Stamford, CT: Thomson Learning, Inc., 2001.

24.

Tan

Zhong

. Adaptive selection of parameters for precise computation of matrix exponential based on Padé approximation. Chinese J Theor Appl Mech 2009; 41: 961–966.

25.

Chapra

Canale

. Numerical methods for engineers. 6th ed. New York: McGraw-Hill, 2009.