Solving differential-algebraic equations by Taylor series (II): Computing the System Jacobian

The authors have developed a Taylor series method for solving numerically an initial-value problem differential-algebraic equation (DAE) that can be of high index, high order, nonlinear, and fully implicit, BIT, 45 (2005), pp. 561–592. Numerical results have shown that this method is efficient and very accurate. Moreover, it is particularly suitable for problems that are of too high an index for present DAE solvers. This paper develops an effective method for computing a DAE’s System Jacobian, which is needed in the structural analysis of the DAE and computation of Taylor coefficients. Our method involves preprocessing of the DAE and code generation employing automatic differentiation. Theory and algorithms for preprocessing and code generation are presented. An operator-overloading approach to computing the System Jacobian is also discussed. AMS subject classification (2000)  34A09, 65L80, 65L05, 41A58     

Breaking the limits: The Taylor series method

This paper discusses several examples of ordinary differential equation (ODE) applications that are difficult to solve numerically using conventional techniques, but which can be solved successfully using the Taylor series method. These results are hard to obtain using other methods such as Runge-Kutta or similar schemes; indeed, in some cases these other schemes are not able to solve such systems at all. In particular, we explore the use of the high-precision arithmetic in the Taylor series method for numerically integrating ODEs. We show how to compute the partial derivatives, how to propagate sets of initial conditions, and, finally, how to achieve the Brouwer’s Law limit in the propagation of errors in long-time simulations. The TIDES software that we use for this work is freely available from a website.     

Numerical solutions of index-1 differential algebraic equations can be computed in polynomial time

The cost of solving an initial value problem for index-1 differential algebraic equations to accuracy ɛ is polynomial in ln(1/ɛ). This cost is obtained for an algorithm based on the Taylor series method for solving differential algebraic equations developed by Pryce. This result extends a recent result by Corless for solutions of ordinary differential equations. The results of the standard theory of information-based complexity give exponential cost for solving ordinary differential equations, being based on a different model.     

SEM Modeling with Singular Moment Matrices Part I: ML-Estimation of Time Series

A structural equation model (SEM) with deterministic intercepts is introduced. The Gaussian likelihood function does not contain determinants of sample moment matrices and is thus well-defined for only one statistical unit. The SEM is applied to the dynamic state space model and compared with the Kalman filter (KF) approach. The likelihood of both methods are shown to be equivalent, but for long time series numerical problems occur in the SEM approach, which are traced to the inversion of the latent state covariance matrix. Both approaches are compared on several aspects. The SEM approach is now open for idiographic (N = 1) analysis and estimation of panel data with correlated units.     

Polynomial cost for solving IVP for high-index DAE

We show that the cost of solving initial value problems for high-index differential algebraic equations is polynomial in the number of digits of accuracy requested. The algorithm analyzed is built on a Taylor series method developed by Pryce for solving a general class of differential algebraic equations. The problem may be fully implicit, of arbitrarily high fixed index and contain derivatives of any order. We give estimates of the residual which are needed to design practical error control algorithms for differential algebraic equations. We show that adaptive meshes are always more efficient than non-adaptive meshes. Finally, we construct sufficiently smooth interpolants of the discrete solution. AMS subject classification (2000) 34A09, 65L80, 68Q25     

Statistical Computations on Biological Rhythms I: Dissecting Variable Cycles and Computing Signature Phases in Activity-Event Time Series

We propose a computational methodology to compute and extract circadian rhythmic patterns from an individual animal’s activity-event time series. This lengthy dataset, composed of a sequential event history, contains an unknown number of latent rhythmic cycles of varying duration and missing waveform information. Our computations aim at identifying the onset signature phase which individually indicates a sharp event intensity surge, where a subject-night ends and a brand new cycle’s subject-day begins, and collectively induces a linearity manifesting the individual circadian rhythmicity and information about the average period. Based on the induced linearity, the least squares criterion is employed to choose an optimal sequence of computed onset signature phases among a finite collection derived from the hierarchical factor segmentation (HFS) algorithm. The multiple levels of coding schemes in the HFS algorithm are designed to extract contrasting patterns of aggregation against sparsity of activity events along the entire temporal axis. This optimal sequence dissects the whole time series into a sequence of rhythmic cycles without model assumptions or ad hoc behavioral definitions regarding the missing waveform information. The performance of our methodology is favorably compared with two popular approaches based on the periodogram in a simulation study and in real data analyses. The computer code and data used in this article are available on the JCGS webpage.     

解一般矩阵方程的GPBiCG(m,l)法

The generalized product bi-conjugate gradient(GPBiCG(m,l))method has been recently proposed as a hybrid variant of the GPBi CG and the Bi CGSTAB methods to solve the linear system Ax=b with non-symmetric coefficient matrix,and its attractive convergence behavior has been authenticated in many numerical experiments.By means of the Kronecker product and the vectorization operator,this paper aims to develop the GPBi CG(m,l)method to solve the general matrix equation■ and the general discrete-time periodic matrix equations■ which include the well-known Lyapunov,Stein,and Sylvester matrix equations that arise in a wide variety of applications in engineering,communications and scientific computations.The accuracy and efficiency of the extended GPBi CG(m,l)method assessed against some existing iterative methods are illustrated by several numerical experiments.     

SEM Modeling with Singular Moment Matrices Part II: ML-Estimation of Sampled Stochastic Differential Equations

Linear stochastic differential equations are expressed as an exact discrete model (EDM) and estimated with structural equation models (SEMs) and the Kalman filter (KF) algorithm. The oversampling approach is introduced in order to formulate the EDM on a time grid which is finer than the sampling intervals. This leads to a simple computation of the nonlinear parameter functionals of the EDM. For small discretization intervals, the functionals can be linearized, and standard software permitting only linear parameter restrictions can be used. However, in this case the SEM approach must handle large matrices leading to degraded performance and possible numerical problems. The methods are compared using coupled linear random oscillators with time-varying parameters and irregular sampling times.     

Optimal solutions of a (3 + 1)-dimensional B-Kadomtsev-Petviashvii equation

We explored and specialized new Lie infinitesimals for the (3 + 1)-dimensional B-Kadomtsev-Petviashvii (BKP) using the commutation product, which results a system of nonlinear ODEs manually solved. Through two stages of Lie symmetry reduction, (3 + 1)-dimensional BKP equation is reduced to nonsolvable nonlinear ODEs using various combinations of optimal Lie vectors. Using the integration and Riccati equation methods, we investigate new analytical solutions for these ODEs. Back substituting to the original variables generates new solutions for BKP. Some selected solutions illustrated through three-dimensional plots.     

Benchmarking with quality-adjusted DEA (Q-DEA) to seek lower-cost high-quality service: Evidence from a U.S.bank application

Benchmarking is a widely cited method to identify and adopt best-practices as a means to improve performance. Data envelopment analysis (DEA) has been demonstrated to be a powerful benchmarking methodology for situations where multiple inputs and outputs need to be assessed to identify best-practices and improve productivity in organizations. Most DEA benchmarking studies have excluded quality, even in service-sector applications such as health care where quality is a key element of performance. This limits the practical value of DEA in organizations where maintaining and improving service quality is critical to achieving performance objectives. In this paper, alternative methods incorporating quality in DEA benchmarking are demonstrated and evaluated. It is shown that simply treating the quality measures as DEA outputs does not help in discriminating the performance. Thus, the current study presents a new, more sensitive, quality-adjusted DEA (Q-DEA), which effectively deals with quality measures in benchmarking. We report the results of applying Q-DEA to a U.S. bank's 200-branch network that required a method for benchmarking to help manage operating costs and service quality. Q-DEA findings helped the bank achieve cost savings and improved operations while preserving service quality, a dimension critical to its mission. New insights about ways to improve branch operations based on the best-practice (high-quality low-cost) benchmarks identified with Q-DEA are also described in the paper. This demonstrates the practical need and potential benefits of Q-DEA and its efficacy in one application, and also suggests the need for further research on measuring and incorporating quality into DEA benchmarking. The review process of this paper was handled by the Edit-in-Chief Peter Hammer.     

Computing eigenspaces with specified eigenvalues of a regular matrix pair (A,B) and condition estimation: theory,algorithms and software

Theory, algorithms and LAPACK-style software for computing a pair of deflating subspaces with specified eigenvalues of a regular matrix pair (A, B) and error bounds for computed quantities (eigenvalues and eigenspaces) are presented. Thereordering of specified eigenvalues is performed with a direct orthogonal transformation method with guaranteed numerical stability. Each swap of two adjacent diagonal blocks in the real generalized Schur form, where at least one of them corresponds to a complex conjugate pair of eigenvalues, involves solving a generalized Sylvester equation and the construction of two orthogonal transformation matrices from certain eigenspaces associated with the diagonal blocks. The swapping of two 1×1 blocks is performed using orthogonal (unitary) Givens rotations. Theerror bounds are based on estimates of condition numbers for eigenvalues and eigenspaces. The software computes reciprocal values of a condition number for an individual eigenvalue (or a cluster of eigenvalues), a condition number for an eigenvector (or eigenspace), and spectral projectors onto a selected cluster. By computing reciprocal values we avoid overflow. Changes in eigenvectors and eigenspaces are measured by their change in angle. The condition numbers yield bothasymptotic andglobal error bounds. The asymptotic bounds are only accurate for small perturbations (E, F) of (A, B), while the global bounds work for all (E, F.) up to a certain bound, whose size is determined by the conditioning of the problem. It is also shown how these upper bounds can be estimated. Fortran 77software that implements our algorithms for reordering eigenvalues, computing (left and right) deflating subspaces with specified eigenvalues and condition number estimation are presented. Computational experiments that illustrate the accuracy, efficiency and reliability of our software are also described.