Adaptive partitioning techniques for ordinary differential equations

Many stiff systems of ordinary differential equations (ODEs) modeling practical problems can be partitioned into loosely coupled subsystems. In this paper the objective of the partitioning is to permit the numerical integration of one time step to be performed as the solution of a sequence of small subproblems. This reduces the computational complexity compared to solving one large system and permits efficient parallel execution under appropriate conditions. The subsystems are integrated using methods based on low order backward differentiation formulas.This paper presents an adaptive partitioning algorithm based on a classical graph algorithm and techniques for the efficient evaluation of the error introduced by the partitioning.The power of the adaptive partitioning algorithm is demonstrated by a real world example, a variable step-size integration algorithm which solves a system of ODEs originating from chemical reaction kinetics. The computational savings are substantial. In memory of Germund Dahlquist (1925–2005).AMS subject classification (2000) 65L06, 65Y05     

One-step 9-stage Hermite–Birkhoff–Taylor DAE solver of order 10

The HBT(10)9 method for ODEs is expanded into HBT(10)9DAE for solving nonstiff and moderately stiff systems of fully implicit differential algebraic equations (DAEs) of arbitrarily high fixed index. A scheme to generate first-order derivatives at off-step points is combined with Pryce scheme which generates high order derivatives at step points. The stepsize is controlled by a local error estimator. HBT(10)9DAE uses only the first four derivatives of y instead of the first 10 required by Taylor’s series method T10DAE of order 10. Dormand–Prince’s DP(8,7)13M for ODEs is extended to DP(8,7)DAE for DAEs. HBT(10)9DAE wins over DP(8,7)DAE on several test problems on the basis of CPU time as a function of relative error at the end of the interval of integration. An index-5 problem is equally well solved by HBT(10)9DAE and T10DAE. On this problem, the error in the solution by DP(8,7)DAE increases as time increases.     

A parareal approach of semi‐linear parabolic equations based on general waveform relaxation

We present a parareal approach of semi‐linear parabolic equations based on general waveform relaxation (WR) at the partial differential equation (PDE) level. An algorithm for initial‐boundary value problem and two algorithms for time‐periodic boundary value problem are constructed. The convergence analysis of three algorithms are provided. The results show that the algorithm for initial‐boundary value problem is superlinearly convergent while both algorithms for the time‐periodic boundary value problem linearly converge to the exact solutions at most. Numerical experiments show that the parareal algorithms based on general WR at the PDE level, compared with the parareal algorithm based on the classical WR at the ordinary differential equations (ODEs) level (the PDEs is discretized into ODEs), require much fewer number of iterations to converge.     

Index minimization of differential-algebraic equations in hybrid analysis for circuit simulation

Modern modeling approaches for circuit analysis lead to differential-algebraic equations (DAEs). The index of a DAE is a measure of the degree of numerical difficulty. In general, the higher the index is, the more difficult it is to solve the DAE. The index of the DAE arising from the modified nodal analysis (MNA) is determined uniquely by the structure of the circuit. Instead, we consider a broader class of analysis method called the hybrid analysis. For linear time-invariant electric circuits, we devise a combinatorial algorithm for finding an optimal hybrid analysis in which the index of the DAE to be solved attains the minimum. The optimal hybrid analysis often results in a DAE with lower index than MNA.      

A verified inexact implicit Runge–Kutta method for nonsmooth ODEs

Structural pounding and oscillations have been extensively investigated by using ordinary differential equations (ODEs). In many applications, force functions are defined by piecewise continuously differentiable functions and the ODEs are nonsmooth. Implicit Runge–Kutta (IRK) methods for solving the nonsmooth ODEs are numerically stable, but involve systems of nonsmooth equations that cannot be solved exactly in practice. In this paper, we propose a verified inexact IRK method for nonsmooth ODEs which gives a global error bound for the inexact solution. We use the slanting Newton method to solve the systems of nonsmooth equations, and interval method to compute the set of matrices of slopes for the enclosure of solution of the systems. Numerical experiments show that the algorithm is efficient for verification of solution of systems of nonsmooth equations in the inexact IRK method. We report numerical results of nonsmooth ODEs arising from simulation of the collapse of the Tacoma Narrows suspension bridge, steel to steel impact experiment, and pounding between two adjacent structures in 27 ground motion records for 12 different earthquakes. This work is partly supported by a Grant-in-Aid from Japan Society for the Promotion of Science and a scholarship from Egyptian Government.     

An analysis of the Prothero–Robinson example for constructing new adaptive ESDIRK methods of order 3 and 4

Explicit singly-diagonally-implicit (ESDIRK) Runge–Kutta methods have usually order reduction if they are applied on stiff ODEs, such as the example of Prothero and Robinson. It can be observed that the numerical order of convergence decreases to the stage order, which is limited to two. In this paper we analyse the Prothero–Robinson example and derive new order conditions to avoid order reduction. New third and fourth order ESDIRK methods are created, which are applied to the Prothero–Robinson example and to an index-2 DAE. Numerical examples show that the new methods have better convergence properties than usual ESDIRK methods.     

Numerical integration of positive linear differential-algebraic systems

In the simulation of dynamical processes in economy, social sciences, biology or chemistry, the analyzed values often represent non-negative quantities like the amount of goods or individuals or the density of a chemical or biological species. Such systems are typically described by positive ordinary differential equations (ODEs) that have a non-negative solution for every non-negative initial value. Besides positivity, these processes often are subject to algebraic constraints that result from conservation laws, limitation of resources, or balance conditions and thus the models are differential-algebraic equations (DAEs). In this work, we present conditions under which both these properties, the positivity as well as the algebraic constraints, are preserved in the numerical simulation by Runge–Kutta or multistep discretization methods. Using a decomposition approach, we separate the dynamic and the algebraic equations of a given linear, positive DAE to give positivity preserving conditions for each part separately. For the dynamic part, we generalize the results for positive ODEs to DAEs using the solution representation via Drazin inverses. For the algebraic part, we use the consistency conditions of the discretization method to derive conditions under which this part of the approximation overestimates the exact solution and thus is non-negative. We analyze these conditions for some common Runge–Kutta and multistep methods and observe that for index-1 systems and stiffly accurate Runge–Kutta methods, positivity is conditionally preserved under similar conditions as for ODEs. For higher index problems, however, none of the common methods is suitable.     

Dynamic iteration may fail for partitioned network simulation

In many technical applications, complex systems are composed of subsystems that are modelled independently using a network approach. This co‐modelling ansatz is completed very naturally by a co‐simulation approach based on dynamic iteration, where the respective DAE submodels are solved iteratively. For partitioned electric network simulation we show that this approach may fail. As alternative we discuss overlapping techniques based on topological network information.     

Differential Elimination–Completion Algorithms for DAE and PDAE

Differential–algebraic equations (DAE) and partial differential–algebraic equations (PDAE) are systems of ordinary equations and PDAEs with constraints. They occur frequently in such applications as constrained multibody mechanics, spacecraft control, and incompressible fluid dynamics.
A DAE has differential index r if a minimum of r +1 differentiations of it are required before no new constraints are obtained. Although DAE of low differential index (0 or 1) are generally easier to solve numerically, higher index DAE present severe difficulties.
Reich et al. have presented a geometric theory and an algorithm for reducing DAE of high differential index to DAE of low differential index. Rabier and Rheinboldt also provided an existence and uniqueness theorem for DAE of low differential index. We show that for analytic autonomous first-order DAE, this algorithm is equivalent to the Cartan–Kuranishi algorithm for completing a system of differential equations to involutive form. The Cartan–Kuranishi algorithm has the advantage that it also applies to PDAE and delivers an existence and uniqueness theorem for systems in involutive form. We present an effective algorithm for computing the differential index of polynomially nonlinear DAE. A framework for the algorithmic analysis of perturbed systems of PDAE is introduced and related to the perturbation index of DAE. Examples including singular solutions, the Pendulum, and the Navier–Stokes equations are given. Discussion of computer algebra implementations is also provided.     

New Algorithm for Computing LQ Suboptimal Output Feedback Gains of Decentralized Control Systems

In this paper, the problem of computing the suboptimal output feedback gains of decentralized control systems is investigated. First, the problem is formulated. Then, the gradient matrices based on the index function are derived and a new algorithm is established based on some nice properties. This algorithm shows that a suboptimal gain can be computed by solving several ordinary differential equations (ODEs). In order to find an initial condition for the ODEs, an algorithm for finding a stabilizing output feedback gain is exploited, and the convergence of this algorithm is discussed. Finally, an example is given to illustrate the proposed algorithm.     

Element-based index reduction in electrical circuit simulation

The numerical simulation of very large scale integrated circuits is an important tool in the development of new industrial circuits. In the course of the last years, this topic has received increasing attention. Common modeling approaches for circuits lead to differential-algebraic systems (DAEs). In circuit simulation, these DAEs are known to have index 2, given some topological properties of the network. This higher index leads to several undesirable effects in the numerical solution of the DAEs. Recent approaches try to lower the index of DAEs to improve the numerical behaviour. These methods usually involve costly algebraic transformations of the equations. Especially, for large scale circuit equations, these transformations become too costly to be efficient. We will present methods that change the topology of the network itself, while replacing certain elements in oder to obtain a network that leads to a DAE of index 1, while not altering the analytical solution of the DAE. This procedure can be performed prior to the actual numerical simulation. The decreasing of the index usually leads to significantly improved numerical behaviour. (© 2006 WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim)     

The Index of an Infinite Dimensional Implicit System

The idea of the index of a differential algebraic equation (DAE) (or implicit differential equation) has played a fundamental role in both the analysis of DAEs and the development of numerical algorithms for DAEs. DAEs frequently arise as partial discretizations of partial differential equations (PDEs). In order to relate properties of the PDE to those of the resulting DAE it is necessary to have a concept of the index of a possibly constrained PDE. Using the finite dimensional theory as motivation, this paper will examine what one appropriate analogue is for infinite dimensional systems. A general definition approach will be given motivated by the desire to consider numerical methods. Specific examples illustrating several kinds of behavior will be considered in some detail. It is seen that our definition differs from purely algebraic definitions. Numerical solutions, and simulation difficulties, can be misinterpreted if this index information is missing.     

WAVEFORM RELAXATION METHODS AND ACCURACY INCREASE

ＷＡＶＥＦＯＲＭ　ＲＥＬＡＸＡＴＩＯＮ　ＭＥＴＨＯＤＳ　ＡＮＤ　ＡＣＣＵＲＡＣＹ　ＩＮＣＲＥＡＳＥＷＡＶＥＦＯＲＭＲＥＬＡＸＡＴＩＯＮＭＥＴＨＯＤＳＡＮＤＡＣＣＵＲＡＣＹＩＮＣＲＥＡＳＥ￥ＳｏｎｇＹｏｎｇｚｈｏｎｇ（ＮａｎｊｉｎｇＮｏｒｍａｌＵｎｉｖ...     

Strict Chebyshev approximation for general systems of linear equations

Summary An ascent exchange algorithm for computing the strict Chebyshev solution to general systems of linear equations is presented. It uses generalized exchange rules to ensure convergence and splits up the entire system into subsystems by means of a canonical decomposition of a matrix obtained by Gaussian elimination methods. All updating procedures are developed and several numerical examples illustrate the efficiency of the algorithm.     

A Structure Preserving Index Reduction Method For MNA

Transient analysis in industrial chip design leads to very large systems of differential-algebraic equations (DAEs). The numerical solution of these DAEs strongly depends on the so called index of the DAE. In general, the higher the index of the DAE is, the more sensitive the numerical solution will be to errors in the computation. So, it is advisable to use mathematical models with small index or to reduce the index. This paper presents an index reduction method that uses information based on the topology of the circuit. In addition, we show that the presented method retains structural properties of the DAE. (© 2006 WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim)     

Order Reduction of General Nonlinear DAE Systems by Automatic Tearing

The object-oriented approach to modelling has recently made possible to build models of large-scale real systems. However, the resulting system of equations is generally a nonlinear DAE (Differential Algebraic Equations) system of large dimension, which must be reduced in some way to make it tractable for numerical solutions. A way to do this is automatic symbolic tearing, which aims at splitting the DAE system into two parts: a core consisting of a reduced implicit DAE system and a set of explicit assignments. The problem is here dealt with by a graph theoretic approach, first proving the NP-completeness in the general case, then formulating the problem with reference to a bipartite graph and finally defining an efficient and flexible algorithm for automatic tearing. It is also shown how the proposed algorithm can easily incorporate both general and domain-specific heuristic rules, and can also be used to deal with equations in vector form. The application to serial multibody systems is considered as a significant example.     

Pryce pre-analysis adapted to some DAE solvers

The ODE solver HBT(12)5 of order 12 (T. Nguyen-Ba, H. Hao, H. Yagoub, R. Vaillancourt, One-step 5-stage Hermite-Birkho-Taylor ODE solver of order 12, Appl. Math. Comput. 211 (2009) 313-328. doi:10.1016/j.amc.2009.01.043), which combines a Taylor series method of order 9 with a Runge-Kutta method of order 4, is expanded into the DAE solver HBT(12)5DAE of order 12. Dormand-Prince’s DP(8, 7)13M is also expanded into the DAE solver DP(8, 7)DAE. Pryce structural pre-analysis, extended ODEs and ODE first-order forms are adapted to these DAE solvers with a stepsize control based on local error estimators and a modified Pryce algorithm to advance integration. HBT(12)5DAE uses only the first nine derivatives of the unknown variables as opposed to the first 12 derivatives used by the Taylor series method T12DAE of order 12. Numerical results show the advantage of HBT(12)5DAE over T12DAE, DP(8, 7)DAE and other known DAE solvers.     

Implicit algorithms for solving the Cauchy problem for the equations of the dynamics of mechanical systems

It is shown that in the numerical solution of the Cauchy problem for systems of second-order ordinary differential equations, when solved for the highest-order derivative, it is possible to construct simple and economical implicit computational algorithms for step-by-step integration without using laborious iterative procedures based on processes of the Newton-Raphson iterative type. The initial problem must first be transformed to a new argument — the length of its integral curve. Such a transformation is carried out using an equation relating the initial parameter of the problem to the length of the integral curve. The linear acceleration method is used as an example to demonstrate the procedure of constructing an implicit algorithm using simple iterations for the numerical solution of the transformed Cauchy problem. Propositions concerning the computational properties of the iterative process are formulated and proved. Explicit estimates are given for an integration stepsize that guarantees the convergence of the simple iterations. The efficacy of the proposed procedure is demonstrated by the numerical solution of three problems. A comparative analysis is carried out of the numerical solutions obtained with and without parametrization of the initial problems in these three settings. As a qualitative test the problem of the celestial mechanics of the “Pleiades” is considered. The second example is devoted to modelling the non-linear dynamics of an elastic flexible rod fixed at one end as a cantilever and coiled in its initial (static) state into a ring by a bending moment. The third example demonstrates the numerical solution of the problem of the “unfolding” of a mechanical system consisting of three flexible rods with given control input.