Linearly implicit quantization-based integration methods for stiff ordinary differential equations

In this paper, new integration methods for stiff ordinary differential equations (ODEs) are developed. Following the idea of quantization-based integration (QBI), i.e., replacing the time discretization by state quantization, the proposed algorithms generalize the idea of linearly implicit algorithms. Also, the implementation of the new algorithms in a DEVS simulation tool is discussed. The efficiency of these new methods is verified by comparing their performance in the simulation of two benchmark problems with that of other numerical stiff ODE solvers. In particular, the advantages of these new algorithms for the simulation of electronic circuits are demonstrated.     

Efficient integration over discontinuities in ordinary differential equation simulations

This paper introduces a new method of detecting and handling discontinuities in arbitrary functions which form part of an ordinary differential equation set. The method has been implemented in conjunction with the Gear[6] integration algorithm for stiff equation sets, published by Hindmarsh[9], but the philosophy applies to any predictor corrector integration algorithm with Nordsieck[10] step size control.     

A type-insensitive ode code based on second derivative formulas

Several effective codes for the solution of stiff ordinary differential equations (ODEs) are based on second derivative formulas. They are inefficient for non-stiff problems. It is shown how to modify such codes to make them reasonably efficient. The modifications make them more efficient for stiff problems as well. Users do not have the information to determine stiffness reliably. The modified codes recognize the type automatically at each step and respond appropriately. It is a great convenience for a user to have one code efficient regardless of the type.     

线性多步法模糊逻辑系统

根据求解常微分方程的线性多步逼近方法构造了线性多步法模糊逻辑系统以辨识由常微分方程所描述的未知动态系统,这种新的模糊逻辑系统在预测状态的同时也能够逼近系统微分方程中的未知函数,文中从理论上对精度作了分析,并用仿真验证了所得的结果。     

Stochastic projective methods for simulating stiff chemical reacting systems

In this paper, stochastic projective methods are proposed to improve the stability and efficiency in simulating stiff chemical reacting systems. The efficiency of existing explicit tau-leaping methods can often severely be limited by the stiffness in the system, forcing the use of small time steps to maintain stability. The methods presented in this paper, namely stochastic projective (SP) and telescopic stochastic projective (TSP) method, can be considered as more general stochastic versions of the recently developed stable projective numerical integration methods for deterministic ordinary differential equations. SP and TSP method are developed by fully re-interpreting and extending the key projective integration steps in the deterministic regime under a stochastic context. These new stochastic methods not only automatically reduce to the original deterministic stable methods when applied to simulating ordinary differential equations, but also carry the enhanced stability property over to the stochastic regime. In some sense, the proposed methods are stochastic generalizations to their deterministic counterparts. As such, SP and TSP method can adopt a much larger effective time step than is allowed for explicit tau-leaping, leading to noticeable runtime speedup. The explicit nature of the proposed stochastic simulation methods relaxes the need for solving any coupled nonlinear systems of equations at each leaping step, making them more efficient than the implicit tau-leaping method with similar stability characteristics. The efficiency benefits of SP and TSP method over the implicit tau-leaping is expected to grow even more significantly for large complex stiff chemical systems involving hundreds of active species and beyond.     

Numerical integration of semidiscrete evolution systems

Methods and algorithms for integrating initial value systems are examined. Of particular interest is efficient and accurate numerical integration of systems of ordinary differential equations that arise on semidiscrete spatial differencing or finite element projection for evolution problems characterized by partial differential equations. Integration schemes for general systems are described. Stiff and oscillatory systems are considered and these motivate selection of specific types of algorithms for certain problem classes. For example, we show that Runge-Kutta methods with extended regions of stability are particularly efficient for moderately stiff dissipative systems derived from parabolic transport equations. The theoretical developments of an earlier paper [1] determine bounds on stiffness and stability and may be used to examine the stiff dissipative or oscillatory nature of the system qualitatively. Order control and stepsize adjustment in variable-order, variable-step algorithms are compared for several integrators applied to stiff and nonstiff initial-value systems arising from representative parabolic evolution problems.     

Linearly implicit time integration methods in real-time applications: DAEs and stiff ODEs

The methods for the dynamical simulation of multi-body systems in real-time applications have to guarantee that the time integration of the equations of motion is always successfully completed within an a priori fixed sampling time interval, typically in the range of 1.0–10.0 ms. Model structure, model complexity and numerical solution methods have to be adapted to the needs of real-time simulation. Standard solvers for stiff and for constrained mechanical systems are implicit and cannot be used straightforwardly in real-time applications because of their iterative strategies to solve the nonlinear corrector equations and because of adaptive strategies for stepsize and order selection. As an alternative, we consider in the present paper noniterative fixed stepsize time integration methods for stiff ordinary differential equations (ODEs) resulting from tree-structured multi-body system models and for differential algebraic equations (DAEs) that result from multi-body system models with loop-closing constraints.     

A New Class of Highly Accurate Solvers for Ordinary Differential Equations

We introduce a new class of predictor-corrector schemes for the numerical solution of the Cauchy problem for non-stiff ordinary differential equations (ODEs), obtained via the decomposition of the solutions into combinations of appropriately chosen exponentials; historically, such techniques have been known as exponentially fitted methods. The proposed algorithms differ from the classical ones both in the selection of exponentials and in the design of the quadrature formulae used by the predictor-corrector process. The resulting schemes have the advantage of significantly faster convergence, given fixed lengths of predictor and corrector vectors. The performance of the approach is illustrated via a number of numerical examples. This work was partially supported by the US Department of Defense under ONR Grant #N00014-07-1-0711 and AFOSR Grants #FA9550-06-1-0197 and #FA9550-06-1-0239.     

Order,stepsize and stiffness switching

The question discussed in this paper is “When should a code switch between stiff and non-stiff options, when should it switch between one order and another and how should it adjust its stepsize from one step to the next.” Criteria are proposed for switching between the options that become available as the integration progresses and these are presented in algorithmic form.     

Robust real-time identification of linear systems with correlated noise

The problem of recursive robust identification of linear discrete-time single-input single-output dynamic systems with correlated disturbances is considered. Problems related to the construction of optimal robust stochastic approximation algorithms in the min-max sense are demonstrated. Since the optimal solution cannot be achieved in practice, several robustified stochastic approximation algorithms are derived on the basis of a suitable non-linear transformation of normalized residuals, as well as step-by-step optimization with respect to the weighting matrix of the algorithm. The convergence of the developed algorithms is established theoretically using the ordinary differential equation approach. Monte Carlo simulation results are presented for the quantitative performance evaluation of the proposed algorithms. The results indicate the most suitable algorithms for applications in engineering practice.     

Nonlinearity problem in the numerical solution of superstiff Cauchy problems

For the numerical solution of Cauchy stiff initial problems, many schemes have been proposed for ordinary differential equation systems. They work well on linear and weakly nonlinear problems. The article presents a study of a number of well-known schemes on nonlinear problems (which include, for example, the problem of chemical kinetics). It is shown that on these problems, the known numerical methods are unreliable. They require a sufficient step reducing at some critical moments, and to determine these moments, sufficiently reliable algorithms have not been developed. It is shown that in the choice of time as an argument, the difficulty is associated with the boundary layer. If the length of the integral curve arc is taken as an argument, difficulties are caused by the transition zone between the boundary layer and regular solution.     

On the simulation of nonlinear bidimensional spiking neuron models

Bidimensional spiking models are garnering a lot of attention for their simplicity and their ability to reproduce various spiking patterns of cortical neurons and are used particularly for large network simulations. These models describe the dynamics of the membrane potential by a nonlinear differential equation that blows up in finite time, coupled to a second equation for adaptation. Spikes are emitted when the membrane potential blows up or reaches a cutoff θ. The precise simulation of the spike times and of the adaptation variable is critical, for it governs the spike pattern produced and is hard to compute accurately because of the exploding nature of the system at the spike times. We thoroughly study the precision of fixed time-step integration schemes for this type of model and demonstrate that these methods produce systematic errors that are unbounded, as the cutoff value is increased, in the evaluation of the two crucial quantities: the spike time and the value of the adaptation variable at this time. Precise evaluation of these quantities therefore involves very small time steps and long simulation times. In order to achieve a fixed absolute precision in a reasonable computational time, we propose here a new algorithm to simulate these systems based on a variable integration step method that either integrates the original ordinary differential equation or the equation of the orbits in the phase plane, and compare this algorithm with fixed time-step Euler scheme and other more accurate simulation algorithms.     

Curvature-based grid step selection for stiff Cauchy problems

A new method of automatic step selection is proposed for the numerical integration of the Cauchy problem for ordinary differential equations. The method is based on using the geometrical characteristics (cuvature and slope) of the integral curve. Formulas have been constructed for the curvature of the integral curve for different choices of multidimensional space. In the two-dimensional case, they turn into well-known formulas, but their general multidimensional form is nontrivial. These formulas have a simple form, are convenient for practical use, and are of independent interest for the differential geometry of multidimensional spaces. For the grids constructed by our method, a procedure of step splitting is proposed that allows one to apply Richardson’s method and to calculate posterior asymptotically precise error estimation for the obtained solution (no such estimates have been found for traditional algorithms of automatic step selection). Therefore, the proposed methods demonstrate significantly superior reliability and validity of the results as compared to calculations by conventional algorithms. In the existing automatic procedures for step selection, steps can be unexpectedly reduced by 2–4 orders of magnitude for no apparent reason. This undermines the reliability of the algorithms. The cause of this phenomenon is explained. The proposed methods are especially effective for highly stiff problems, which is illustrated by examples of calculations.     

A computationally efficient technique for transient analysis of repairable Markovian systems

A technique to design efficient methods using a combination of explicit (non-stiff) and implicit (stiff) ODE methods for numerical transient analysis of repairable Markovian systems is proposed. Repairable systems give rise to stiff Markov chains due to extreme disparity between failure rates and repair rates. Our approach is based on the observation that stiff Markov chains are non-stiff for an initial phase of the solution interval. A non-stiff ODE method is used to solve the model for this phase and a stiff ODE method is used to solve the model for the rest of the duration until the end of solution interval. A formal criterion to determine the length of the non-stiff phase is described. A significant outcome of this approach is that the accuracy requirement automatically becomes a part of model stiffness. Two specific methods based on this approach have been implemented. Both the methods use the Runge-Kutta-Fehlberg method as the non-stiff method. One uses the TR-BDF2 method as the stiff method while the other uses an implicit Runge-Kutta method as the stiff method. Numerical results obtained from solving dependability models of a multiprocessor system and an interconnection network are presented. These results show that the methods obtained using this approach are much more efficient than the corresponding stiff methods which have been proposed to solve stiff Markov models.     

Recursive and Residual Algorithms for the Efficient Numerical Integration of Multi-Body Systems

In this paper a family of methods for multi-body dynamic simulation is introduced. Equations of motion are obtained using a set of Cartesian coordinates and projected onto a set of independent relative coordinates using the concept of velocity transformation. Open-chain systems are solved directly following either a fully recursive or a semi-recursive procedure. Closed-chain systems are solved in two steps; kinematic loops are opened by removing either some kinematic joints or a rigid body, and the resulting open-chain system is solved; closure-of-the-loop conditions are imposed by means of a second velocity transformation. The dynamic formalisms have been developed so as to handle both non-stiff and stiff systems. Non-stiff systems are solved by means of an Adams–Bashforth–Moulton numerical integration scheme, which requires the computation of the function derivatives. Stiff problems are integrated by using either BDF or NDF methods, which require the computation of the residual of the equations of motion and, optionally, the evaluation of the Jacobian matrix. The proposed algorithms have been implemented using an Object-Oriented Programming approach that makes it possible to re-use the source code, keeping programs smaller, cleaner and easier to maintain. Practical examples that illustrate the performance of these implementations are included. These examples have also been solved using a commercial multi-body simulation package and comparative results are included. In most cases, the algorithms here presented outperform those implemented in the commercial package, leading to important savings in terms of total computation times.     

Approximate conditional least squares estimation of a nonlinear state-space model via an unscented Kalman filter

The problem of estimating a nonlinear state-space model whose state process is driven by an ordinary differential equation (ODE) or a stochastic differential equation (SDE), with discrete-time data is studied. A new estimation method is proposed based on minimizing the conditional least squares (CLS) with the conditional mean function computed approximately via the unscented Kalman filter (UKF). Conditions are derived for the UKF–CLS estimator to preserve the limiting properties of the exact CLS estimator, namely, consistency and asymptotic normality, under the framework of infill asymptotics, i.e. sampling is increasingly dense over a fixed domain. The efficacy of the proposed method is demonstrated by simulation and a real application.     

Dynamics and control of vapor recompression distillation

Vapor recompression distillation is an energy integrated distillation configuration which works on the principle of a heat pump. The tight material and energy integration in such columns shows a potential for intricate dynamics. In this paper, a systematic modeling framework is developed which explicitly captures the discrepancies between different material and energy flows present in such columns. This feature is documented to lead to stiff dynamic equations and multiple time-scale dynamics. Through a nested application of singular perturbations, reduced order non-stiff models are derived which capture the dynamics in each time-scale. A hierarchical control strategy is then proposed exploiting this time-scale multiplicity. A case of propane–propylene separation is considered to illustrate these results and demonstrate the effectiveness of the proposed control strategy via a simulation case study.     

Solving Initial Value Problems for Ordinary Differential Equations by two approaches: BDF and piecewise-linearized methods

Many scientific and engineering problems are described using Ordinary Differential Equations (ODEs), where the analytic solution is unknown. Much research has been done by the scientific community on developing numerical methods which can provide an approximate solution of the original ODE. In this work, two approaches have been considered based on BDF and Piecewise-linearized Methods. The approach based on BDF methods uses a Chord-Shamanskii iteration for computing the nonlinear system which is obtained when the BDF schema is used. Two approaches based on piecewise-linearized methods have also been considered. These approaches are based on a theorem proved in this paper which allows to compute the approximate solution at each time step by means of a block-oriented method based on diagonal Padé approximations. The difference between these implementations is in using or not using the scale and squaring technique.Five algorithms based on these approaches have been developed. MATLAB and Fortran versions of the above algorithms have been developed, comparing both precision and computational costs. BLAS and LAPACK libraries have been used in Fortran implementations. In order to compare in equality of conditions all implementations, algorithms with fixed step have been considered. Four of the five case studies analyzed come from biology and chemical kinetics stiff problems. Experimental results show the advantages of the proposed algorithms, especially when they are integrating stiff problems.     

带非同位干扰执行动态有限维系统的输出追踪

本文考虑了边界带有非同位干扰的波动方程和常微分方程串联系统的性能输出追踪问题. 边界带有非同位干扰的波动方程可以看作常微分方程控制系统的执行动态. 通过利用轨道规划方法将非同位干扰变换到控制通道, 于是克服了非同位干扰带来的困难. 与现有Backstepping方法不同, 文章利用动态补偿法给出了一个新的状态反馈控制器, 使得控制器设计以及所设计的控制器本身更加简单, 证明了闭环系统解的一致有界性与性能输出的指数追踪. 数值模拟表明, 所给的方法是非常有效的     

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.