Precise time integration for linear two-point boundary value problems

https://doi.org/10.1016/j.amc.2005.08.001Get rights and content

Abstract

The precise time integration (PTI) method is proposed to solve the linear two-point boundary value problem (TPBVP). By employing the method of dimensional expanding, the non-homogeneous ordinary differential equations (ODEs) can be transformed into homogeneous ones and then the original PTI algorithms can be applied directly to the TPBVP. The PTI consists of two methods: the method of matrix exponential and the method of Riccati equations, both of which utilize the merit of 2N algorithm and guarantee high precise numerical results. The method of matrix exponential follows the similar scheme proposed for initial value problem of ODEs, and uses the matrix exponential to link boundary conditions between the two points. The method of Riccati equations employs the Riccati equations to express the relationships of two point boundary conditions. And then in terms of the relationships of boundary conditions at two points, the full conditions at initial point can be obtained and then the TPBVP can be transformed into initial value problem and then solved by direct time marching scheme. With some modifications, the above algorithms can be directly extended to the infinite-interval problem and the variant coefficient ODEs problem. In the program implementation, the object-oriented (OO) design of PTI is proposed to demonstrate the applicability and easy maintenance of OO techniques in numerical computation. Finally, four selected numerical examples are given to show the high precision characteristics of PTI.

Introduction

To solve two point boundary value problems (TPBVP) of ordinary differential equations (ODEs), various methods, such as the initial parameter method and the shooting method, were derived [1]. However these methods are not always effective and may cause ill-conditioning for a problem with long interval. The precise time integration method (PTI) was firstly proposed by Zhong and Williams [2] for linear initial value problem of structural dynamics. Because of its prominent numerical advantages, such as high precision, zero-amplitude rate of decay, zero-period specific elongation and non-overstep, PTI has attracted many attentions during the last decade and its application has been broadened [3], [4], [5], [6], [7], [8], [9]. In addition to initial value problem, Zhong continued to extend the PTI to the TPBVP by using the Riccati equations. In Ref [10], [11], the PTI was applied to solve the TPBVP for homogeneous ODEs. In conjunction with the finite strip method, PTI was employed to solve the TPBVP of elastic plate [12], where the non-homogeneous terms in ODEs are constant. Detailed review on PTI method was given by Zhong [13].

In fact, the non-homogeneous terms in ODEs should not become an obstacle for the application of PTI. In Ref. [14], by using the method of dimensional expanding, the non-homogeneous ODEs were transformed into homogeneous ones and the PTI for initial value problem was directly applied successfully while avoiding the inverse matrix computation. The time-dependent n-dimensional ODEs system are expressed asv˙=Hv+r,where v, H, r are the n-dimensional unknown system state variable, the n × n constant matrix coefficient and the n-dimensional time-dependent non-homogeneous vector. If r is a linear combination of f and f satisfies m-dimensional homogeneous ODEs, then we haver=Hvff,f˙=Hfff,where Hvf and Hff are n × m and m × m matrices, respectively. Then Eq. (1) can be rewritten asv˙=Hv,wherev=vf,H=HHvf0Hff.The idea of dimensional expanding and the transformation from non-homogeneous ODEs to homogeneous ODEs is very simple but the benefits in numerical computation are remarkable.

Here the key issue is how to determine Eq. (2a), (2b). In fact, the expression of f in Eq. (2a), (2b) has a broad range to select. For example, polynomial functions, sinusoidal functions, Fourier series can satisfy Eq. (2b). If r is expressed as a combination of the above functions or series, we can then obtain explicit expression of Eq. (2a), (2b). However, under extreme situation where we cannot determine the explicit expression of Eq. (2a), (2b) and r directly, we can use polynomial series or Fourier series to approximate actual r and then use the series to obtain Eq. (2a), (2b). In general, such an approximation is effective. Here, f is referred to as the basic modes of non-homogeneous vector.

The TPBVP discussed in this paper is formulated asq˙p˙=ADBCqp+rq(t)rp(t),where q, p are vectors with nq and np dimensions (n = nq + np), respectively. Then A, B, C, D, rq(t) and rp(t) are the corresponding matrix/vector according to Eq. (1). The conditions of the ODEs are prescribed at initial time t0 and the final time tf asq(t=t0)=q0,p(t=tf)=pf.The TPBVP with dimensional expanding remains the same as Eq. (3) while the specific expressions of matrix and vector should be determined according to Eq. (5).

In this paper, our focus will be on development of solution processes for Eq. (5). Firstly, the method of matrix exponential and the method of Riccati equations will be discussed in details in Sections 2 Method of matrix exponential, 3 Method of Riccati equations for the ODEs problems with constant coefficients. And then the extension of PTI method to infinite interval problem and variant coefficient problem are presented in Section 4. In the program implementation, the object-oriented design for PTI will be proposed in order to fully exploit the re-usability and easy extensibility of object-oriented techniques in numerical computation. Finally, in Section 6, four selected examples are employed to demonstrate the remarkable advantages of PTI.

Section snippets

Method of matrix exponential

In the conventional PTI, the time marching scheme for initial value problem can be expressed asvk+1=Tvk+0sexp[H(s-s˜)]·r(tk+s˜)ds˜,where s is the time step length and T is the matrix exponential,T=exp(Hs).The key feature of PTI is the precise computation of matrix exponential which was detailed in Ref. [2], [3], [13], [14].

In particular, when r varies linearly with the time interval [tk, tk+1],r(t)=r0+r1(t-tk).Eq. (7) can take the formvk+1=T[vk+H-1(r0+H-1r1)]-H-1[r0+H-1r1+r1s].The expressions

Method of Riccati equations

In Refs. [10], [11], the PTI with Riccati equations was employed to solve the homogeneous ODEs. By using the method of dimensional expanding, Eq. (5) can be rewritten asq˙p˙=ADBCqp,whereq˙=qf,A=AHqvf0Hff,B=BHpvf,D=D0.In order to deduce the relationship connecting q0, p0 at t0 and qf, pf at tf, we first consider the sub-interval [ta, tb], t0  ta < tb  tf. Then the relationship between qa, pa at ta and qb, pb at tb can be formulated asqbpa=F-ΓQEqapb,whereF=FFvf0Fff,Γ=Γ0,Q=QQvf.

Extension of PTI methods

The above methods for constant coefficient problem can be directly extended to infinite interval and variant coefficient problems.

Object-oriented program implementation

During the last decade, the object-oriented (OO) techniques have been applied to numerical computations, in particular the recently developed OO finite element methods [16]. The OO programming provides the concepts of abstract data type, class hierarchies, inheritance and polymorphism for expressing common behavior and operations on data. Therefore the application of OO techniques to numerical computation enables not only to exploit the full re-usability and easy extendibility of the codes, but

Numerical examples

In the following examples, the algorithm parameters of PTI are selected as N = 20 and p = 4.

Example 1

A two dimensional stiff problem, adopted from Ref. [1, p. 739], is formulated asdu/dtdv/dt=9981998-999-1999uv+r(t)r(t).The stiff characteristic is that the great difference of two eigenvalues of the problem λ1 = −1 and λ2 = −1000. In Ref. [1], the problem was expressed as an initial problem with the initial conditions: u(0) = 1 and v(0) = 0. When the f(t) is specified, it is possible to find the exact solution, say u

Conclusions

Two point boundary value problem (TPBVP) is very important in engineering, and in some problem, high precise numerical results are necessary. The precise time integration (PTI) provides a means to achieve this goal. The conclusions of the present research are summarized as follows:

  • (a)

    By using the method of dimensional expanding, two PTI methods, the method of matrix exponential and the method of Riccati method are easily extended to solve TPBVP with non-homogeneous terms.

  • (b)

    By using the 2N algorithm,

Acknowledgement

The authors are grateful to the financial supports of the National Science Foundation of China (NSFC project 10421002), the Australian Research Council (ARC Linkage International Award LX 0348548), and the NSFC-ARC international collaboration project.

References (20)

There are more references available in the full text version of this article.

Cited by (22)

  • Simulating nonlinear aeroelastic responses of an airfoil with freeplay based on precise integration method

    2015, Communications in Nonlinear Science and Numerical Simulation
    Citation Excerpt :

    These models were generally solved by time-marching integration techniques such as the Runge–Kutta (RK) and Newmark methods, etc [8–11]. The precise integration method (PIM) initiated by Zhong [12] has been widely applied to various problems modeled by ordinary differential equations such as structural dynamics, optimal control, and flexible multi-body dynamics problems, etc [13–17]. This method is famous for its high accuracy and computation efficiency.

  • The development and extended applications of precise integration method

    2024, Jisuan Lixue Xuebao/Chinese Journal of Computational Mechanics
  • Novel Global Harmony Search Algorithm for Nonlinear Two-Point Boundary Value Problem

    2021, Proceedings of the 33rd Chinese Control and Decision Conference, CCDC 2021
View all citing articles on Scopus
View full text