Practical aspects of numerical time integration

This paper examines some practical aspects in the numerical time integration of structural equations of motion. The concept of stability region in the complex frequency/time-step (ω · Δt) plane is proposed as the basis of stability evaluation, as opposed to the concept of conditional (or unconditional) stability. Stabilization of numerical computations via the introduction of artificial damping and the composition of existing schemes is discussed. It is shown that stabilization by artificial damping fails for any explicit scheme contrary to the notion that a frequency-proportional damping would suppress the high-frequency components. It is also shown that among explicit schemes examined in this paper, the central difference scheme is preferred from the stability considerations. For nonlinear problems, the linearly extrapolated pseudo-force solution procedure is adopted to assess how nonlinearities affect the stability of implicit schemes. Finally, the impact of stability and accuracy characteristics upon numerical computation is discussed.     

Energy consistent algorithms for dynamic finite deformation plasticity

This paper addresses the theoretical development and numerical implementation of energy consistent algorithms for dynamic elastoplasticity, emphasizing finite strain constitutive formulations so that unconditional stability of the algorithms is assured even in the fully nonlinear regime. The key concept behind energy consistency is the requirement that the discretized system obey an a priori stability estimate, which may be derived in general using the first and second laws of thermodynamics. This approach to computational dynamic plasticity differs from typical application of traditional algorithms (such as Newmark or Hilber–Hughes–Taylor-α methods), where local time integration schemes for plasticity laws are developed somewhat independently from the global time integration scheme for the equations of motion, without explicit consideration of thermodynamical restrictions. Two algorithms based on both additive and multiplicative finite deformation plasticity model are formulated within the energy consistent framework. Both algorithms possess the desirable feature of nonlinear stability of previous energy–momentum algorithms for elastodynamics.     

Computational Schemes for Flexible,Nonlinear Multi-Body Systems

This paper deals with the development of computational schemes for the dynamic analysis of flexible, nonlinear multi-body systems. The focus of the investigation is on the derivation of unconditionally stable time integration schemes for these types of problem. At first, schemes based on Galerkin and time discontinuous Galerkin approximations applied to the equations of motion written in the symmetric hyperbolic form are proposed. Though useful, these schemes require casting the equations of motion in the symmetric hyperbolic form, which is not always possible for multi-body applications. Next, unconditionally stable schemes are proposed that do not rely on the symmetric hyperbolic form. Both energy preserving and energy decaying schemes are derived that both provide unconditionally stable schemes for nonlinear multi-body systems. The formulation of beam and flexible joint elements, as well as of the kinematic constraints associated with universal and revolute joints. An automated time step selection procedure is also developed based on an energy related error measure that provides both local and global error levels. Several examples of simulation of realistic multi-body systems are presented which illustrate the efficiency and accuracy of the proposed schemes, and demonstrate the need for unconditional stability and high frequency numerical dissipation.     

Splitting Methods for Three-Dimensional Transport Models with Interaction Terms

We investigate the use of splitting methods for the numerical integration of three-dimensional transport-chemistry models. In particular, we investigate various possibilities for the time discretization that can take advantage of the parallelization and vectorization facilities offered by multi-processor vector computers. To suppress wiggles in the numerical solution, we use third-order, upwind-biased discretization of the advection terms, resulting in a five-point coupling in each direction. As an alternative to the usual splitting functions, such as co-ordinate splitting or operator splitting, we consider a splitting function that is based on a three-coloured hopscotch-type splitting in the horizontal direction, whereas full coupling is retained in the vertical direction. Advantages of this splitting function are the easy application of domain decomposition techniques and unconditional stability in the vertical, which is an important property for transport in shallow water. The splitting method is obtained by combining the hopscotch-type splitting function with various second-order splitting formulae from the literature. Although some of the resulting methods are highly accurate, their stability behaviour (due to horizontal advection) is quite poor. Therefore we also discuss several new splitting formulae with the aim to improve the stability characteristics. It turns out that this is possible indeed, but the price to pay is a reduction of the accuracy. Therefore, such methods are to be preferred if accuracy is less crucial than stability; such a situation is frequently encountered in solving transport problems. As part of the project TRUST (Transport and Reactions Unified by Splitting Techniques), preliminary versions of the schemes are implemented on the Cray C98 4256 computer and are available for benchmarking.     

Implicit–Explicit Schemes for BGK Kinetic Equations

In this work a new class of numerical methods for the BGK model of kinetic equations is presented. In principle, schemes of any order of accuracy in both space and time can be constructed with this technique. The methods proposed are based on an explicit–implicit time discretization. In particular the convective terms are treated explicitly, while the source terms are implicit. In this fashion even problems with infinite stiffness can be integrated with relatively large time steps. The conservation properties of the schemes are investigated. Numerical results are shown for schemes of order 1, 2 and 5 in space, and up to third-order accurate in time.     

Time Explicit Schemes and Spatial Finite Differences Splittings

In this article, we conjugate time marching schemes with Finite Differences splittings into low and high modes in order to build fully explicit methods with enhanced temporal stability for the numerical solutions of PDEs. The main idea is to apply explicit schemes with less restrictive stability conditions to the linear term of the high modes equation, in order that the allowed time step for the temporal integration is only determined by the low modes. These conjugated schemes were developed in [10] for the spectral case and here we adapt them to the Finite Differences splittings provided by Incremental Unknowns, which steems from the Inertial Manifolds theory. We illustrate their improved capabilities with numerical solutions of Burgers equations, with uniform and nonuniform meshes, in dimensions one and two, when using modified Forward–Euler and Adams–Bashforth schemes. The resulting schemes use time steps of the same order of those used by semi-implicit schemes with comparable accuracy and reduced computational costs.     

A family of stable numerical solvers for the shallow water equations with source terms

In this work we introduce a multiparametric family of stable and accurate numerical schemes for 1D shallow water equations. These schemes are based upon the splitting of the discretization of the source term into centered and decentered parts. These schemes are specifically designed to fulfill the enhanced consistency condition of Bermúdez and Vázquez, necessary to obtain accurate solutions when source terms arise. Our general family of schemes contains as particular cases the extensions already known of Roe and Van Leer schemes, and as new contributions, extensions of Steger–Warming, Vijayasundaram, Lax–Friedrichs and Lax–Wendroff schemes with and without flux-limiters. We include some meaningful numerical tests, which show the good stability and consistency properties of several of the new methods proposed. We also include a linear stability analysis that sets natural sufficient conditions of stability for our general methods.     

On the precise integration methods based on Padé approximations

The precise time step integration method proposed for linear time-invariant homogeneous dynamic systems can give precise numerical results approaching exact solutions at the integration points, but it is conditionally stable. By combining Padé approximation and the generalized Padé approximation of the matrix exponential function in precise integration, and by using three different types of quadrature formulae, a new generalized family of precise time step integration methods is developed to achieve unconditional stability and arbitrary order of accuracy. Numerical studies indicate that they are unconditionally stable algorithms with controllable numerical dissipation. They also demonstrate the validity and efficiency of these algorithms.     

Application of Runge–Kutta–Rosenbrock Methods to the Analysis of Flexible Multibody Systems

Numerical integration methods are discussed for general equations of motion for multibody systems with flexible parts, which are fairly stiff, time-dependent and non-linear. A family of semi-implicit methods, which belong to the class of Runge–Kutta–Rosenbrock methods, with rather weak non-linear stability properties, are developed. These comprise methods of first, second and third order of accuracy that are A-stable and L-stable and hence introduce numerical damping and the filtering of high frequency components. It is shown, both from theory and examples, that it is generally preferable to use deformation mode coordinates to global nodal coordinates as independent variables in the formulation of the equations of motion. The methods are applied to a series of examples consisting of an elastic pendulum, a beam supported by springs, a four-bar mechanism, and a robotic manipulator with collocated control.     

Unconditional Stability, Monotonicity and Accuracy of the Linear Exponential Interpolation Function For the Burgers Equation

The linear exponential function used as an interpolant in the generation of finite difference representations of convection-diffusion equations, is shown to possess several desirable properties not normally found in representations derived from polynomials. These properties include unconditional stability if used in implicit methods, nonoscillatory solutions without the use of flux limiters, and high accuracy. It is also shown that all of these properties are maintained to second-order if the function is used as a spline basis. Comparisons with a number of other schemes are given using the nonlinear Burgers equation. The work here provides the background for the functions arising from a new generalization of Roscoe's method which will be presented in another paper.     

Numerical study of one-step explicit-implicit methods of L-equivalent stiffly accurate two-stage Runge-Kutta schemes

The paper describes one-step methods for numerical integration of the Cauchy problem for systems of ordinary differential equations free from iterations and coinciding on linear problems (autonomous and non-autonomous) with stiffly accurate implicit two-stage Runge-Kutta (RK) schemes. The numerical study of their accuracy is performed on stiff tests, i.e., the autonomous Kaps system and non-autonomous Protero-Robinson problem.     

Numerical integration in the reconstruction of cardiac action potentials using Hodgkin-Huxley-type models

A comparison between traditional numerical integration methods and a new hybrid integration method for the reconstruction of action potential activity is presented, using a mathematical model of the cardiac Purkinje fiber (MNT model). It is shown that the hybrid integration method reduces importantly the overall computation time required for solving the Hodgkin-Huxley differential equations describing membrane electrical events. To accomplish this, the particular form of the gating variable equations is exploited to reformulate the step-by-step computation. In this way, the time increment can be made much larger compared with traditional methods when the membrane potential changes slowly. A mathematical analysis of the hybrid integration method is presented also, together with a numerical verification of its performance both for the propagated and nonpropagated membrane action potential. It is shown that the local error, that is the error arising at each integration step, and the cumulative integration error are strictly controlled by the membrane potential offset. Using the MNT model, the nonpropagated cardiac Purkinje action potential can be reconstructed in real time with an accuracy of 1% for the potential and 5% for the time of occurrence of its main features. In reconstructing propagated events, the hybrid integration method allows computation time savings by a factor of 10 or more compared to accurate Runge-Kutta schemes.     

On the convergence of compact difference schemes

Difference schemes that are compact in space, i.e., schemes constructed on a two- or three-point stencil in each spatial direction, are more efficient and convenient for boundary condition formulation than other high-order accurate schemes. Originally, these schemes were developed primarily to obtain smooth solutions. In the last two decades, compact schemes have been actively used to compute gas dynamic flows with shock waves. However, when a numerical solution with guaranteed accuracy is desired, the actual properties of difference schemes have to be known in the calculation of solutions with discontinuities. For some widely used compact schemes, this issue has not yet been well studied. The properties of compact schemes constructed by the method of lines are examined in this paper. An initial-boundary value problem for the linear heat equation with discontinuous initial data is used as a test problem. In the method of lines, the spatial derivative in the heat equation is approximated on a two-point stencil according to a fourth-order accurate compact differentiation formula. The resulting evolution system of ordinary differential equations is solved using various implicit one-step two- and three-stage schemes of the second and third order of accuracy. The relation between the properties of the stability function of a scheme and the spatial monotonicity of the numerical solution is analyzed. In computations over long time intervals, the compact schemes are shown to be superior to traditional schemes based on the second-order accurate three-point approximation of the spatial derivative.     

Concepts and Application of Time-Limiters to High Resolution Schemes

A new class of implicit high-order non-oscillatory time integration schemes is introduced in a method-of-lines framework. These schemes can be used in conjunction with an appropriate spatial discretization scheme for the numerical solution of time dependent conservation equations. The main concept behind these schemes is that the order of accuracy in time is dropped locally in regions where the time evolution of the solution is not smooth. By doing this, an attempt is made at locally satisfying monotonicity conditions, while maintaining a high order of accuracy in most of the solution domain. When a linear high order time integration scheme is used along with a high order spatial discretization, enforcement of monotonicity imposes severe time-step restrictions. We propose to apply limiters to these time-integration schemes, thus making them non-linear. When these new schemes are used with high order spatial discretizations, solutions remain non-oscillatory for much larger time-steps as compared to linear time integration schemes. Numerical results obtained on scalar conservation equations and systems of conservation equations are highly promising.     

Numerical Approximations for the Cahn–Hilliard Phase Field Model of the Binary Fluid-Surfactant System

In this paper, we consider the numerical approximations for the commonly used binary fluid-surfactant phase field model that consists two nonlinearly coupled Cahn–Hilliard equations. The main challenge in solving the system numerically is how to develop easy-to-implement time stepping schemes while preserving the unconditional energy stability. We solve this issue by developing two linear and decoupled, first order and a second order time-stepping schemes using the so-called “invariant energy quadratization” approach for the double well potentials and a subtle explicit-implicit technique for the nonlinear coupling potential. Moreover, the resulting linear system is well-posed and the linear operator is symmetric positive definite. We rigorously prove the first order scheme is unconditionally energy stable. Various numerical simulations are presented to demonstrate the stability and the accuracy thereafter.     

Algorithmic Foundations of Creation of an Intelligent Software Tool for Investigation and Solution of Cauchy Problems for Systems of Ordinary Differential Equations

Methods are proposed for computer investigation of properties of systems of ordinary differential equations. Based on these methods, algorithms are created for computation of the value of the integration step that provides the stability of a numerical method and obtaining its results with a preassigned accuracy. The components and modes of operation of an intelligent software tool supporting the proposed methods are presented.     

Numerical Methods and Simulations for the Dynamics of One-Dimensional Zakharov–Rubenchik Equations

In this paper, we propose and study several accurate numerical methods for solving the one-dimensional Zakharov–Rubenchik equations (ZRE). We begin with a review on the important properties of the ZRE, including the solitary wave solutions and the various conservation laws. Then we propose a very efficient and accurate numerical method based on the time-splitting technique and the Fourier pseudo-spectral (TSFP) method. Next, we propose some conservative and non-conservative types of finite difference time domain methods, including a Crank–Nicolson finite difference method that conserves the mass and the energy of the system in the discrete level. Discrete conservation laws and numerical stability of all the proposed methods are analyzed. Comparisons between different methods in the efficiency, stability and accuracy are carried out, which identifies that the TSFP method is the most efficient and accurate numerical method among all the methods. Lastly, we apply the TSFP method to simulate and study the dynamics of the solitons in the ZRE numerically.     

Quasi-wavelet method for time-dependent fractional partial differential equation

In this decade, many new applications in engineering and science are governed by a series of fractional partial differential equations. In this paper, we propose a novel numerical method for a class of time-dependent fractional partial differential equations. The time variable is discretized by using the second order backward differentiation formula scheme, and the quasi-wavelet method is used for spatial discretization. The stability and convergence properties related to the time discretization are discussed and theoretically proven. Numerical examples are obtained to investigate the accuracy and efficiency of the proposed method. The comparisons of the present numerical results with the exact analytical solutions show that the quasi-wavelet method has distinctive local property and can achieve accurate results.     

Toward Non Commutative Numerical Analysis: High Order Integration in Time

We review some methods for high precision time integration: it is not easy to ensure stability, precision and numerical efficiency at the same time. Operator splitting—when it works—can be a good way to satisfy all these constraints; in some cases, the order of the splitting schemes can be enhanced by extrapolation; nevertheless, the applicability of splitting is limited due to non commutativity. As an alternative to splitting, we introduce preconditioned Runge–Kutta (PRK) schemes: the preconditioning is included in the scheme, instead of being put aside for implementation. Examples of PRK schemes are given including the extrapolation of the residual smoothing scheme, and sufficient conditions for stability are described.