A Novel Algorithm for the Numerical Integration of Systems of Ordinary Differential Equations Arising in Chemical Problems

Simulation of large networks of chemical reactions via the numerical integration of large systems of ordinary differential equations is of growing importance in real-world problems. We propose an attractive novel numerical integration method, that is largely independent from ill-conditioning and is suitable for any nonlinear problem; moreover, the method, being exact for linear problems, is especially precise for quasi-linear problems, the most frequent kind in the real world. The method is based on a new approach to the computation of a matrix exponential, includes an automatic correction of rounding errors, is not too expensive computationally, and lends itself to a short and robust software implementation that can be easily inserted in large simulation packages. A preliminary numerical verification has been performed, with encouraging results, on two sample problems. The full source listing (in standard C language) of an academic version of the algorithm is freely available on request (e-mail address: Valerio.Parisi@roma2.infn.it), together with a very simple but very stiff chemical problem.     

Wavelet-based collocation method for stiff systems in process engineering

Abrupt phenomena in modelling real-world systems such as chemical processes indicate the importance of investigating stiff systems. However, it is difficult to get the solution of a stiff system analytically or numerically. Two such types of stiff systems describing chemical reactions were modelled in this paper. A numerical method was proposed for solving these stiff systems, which have general nonlinear terms such as exponential function. The technique of dealing with the nonlinearity was based on the Wavelet-Collocation method, which converts differential equations into a set of algebraic equations. Accurate and convergent numerical solutions to the stiff systems were obtained. We also compared the new results to those obtained by the Euler method and 4th order Runge–Kutta method.     

Fast and reliable numerical methods to simulate complex chemical kinetic mechanisms

Models that simulate atmospheric photochemistry require the use of a stiff ordinary differential equations (ODEs) solver. Since the simulation of the chemical transformations taking place in the system takes up to 80 percent of the CPU time, the numerical solver must be computationally fast. Also, the residual error from the solver must be small. Because most accurate solvers are relatively slow, modelers continue to search for timely, yet accurate integration methods. Over the past years an extensive number of articles have been dedicated to this subject. One of the highly debated questions is whether one should construct specialized algorithms or instead use general methods for stiff ODEs. In the present article we use the second alternative. We apply three linearly (semi-)implicit methods from the classical stiff ODE literature which we modified to implement the sparse routines to solve the system of equations describing a complex kinetic mechanism. © 1998 John Wiley & Sons, Inc. Int J Chem Kinet 30: 349–358, 1998     

A combined explicit-implicit method for high accuracy reaction path integration

We present the use of an optimal combined explicit-implicit method for following the reaction path to high accuracy. This is in contrast to most purely implicit reaction path integration algorithms, which are only efficient on stiff ordinary differential equations. The defining equation for the reaction path is considered to be stiff, however, we show here that the reaction path is not uniformly stiff and instead is only stiff near stationary points. The optimal algorithm developed in this work is a combination of explicit and implicit methods with a simple criterion to switch between the two. Using three different chemical reactions, we combine and compare three different integration methods: the implicit trapezoidal method, an explicit stabilized third order algorithm implemented in the code DUMKA3 and the traditional explicit fourth order Runge-Kutta method written in the code RKSUITE. The results for high accuracy show that when the implicit trapezoidal method is combined with either explicit method the number of energy and gradient calculations can potentially be reduced by almost a half compared with integrating either method alone. Finally, to explain the improvements of the combined method we expand on the concepts of stability and stiffness and relate them to the efficiency of integration methods.     

A cutoff phenomenon in accelerated stochastic simulations of chemical kinetics via flow averaging (FLAVOR-SSA)

We present a simple algorithm for the simulation of stiff, discrete-space, continuous-time Markov processes. The algorithm is based on the concept of flow averaging for the integration of stiff ordinary and stochastic differential equations and ultimately leads to a straightforward variation of the the well-known stochastic simulation algorithm (SSA). The speedup that can be achieved by the present algorithm [flow averaging integrator SSA (FLAVOR-SSA)] over the classical SSA comes naturally at the expense of its accuracy. The error of the proposed method exhibits a cutoff phenomenon as a function of its speed-up, allowing for optimal tuning. Two numerical examples from chemical kinetics are provided to illustrate the efficiency of the method.     

Numerical Techniques for the Solution of the Molecular Weight Distribution in Polymerization Mechanisms,State of the Art

The molecular weight distribution (MWD) is possibly the most important characteristic of a polymer. Polymers derive many of their physical properties from their MWD. Therefore, since the origins of polymer science, the theory provides a link between the kinetic mechanism and the mathematical expression of the MWD, and there are analytical solutions for ideal cases. However, the MWD formed in real‐life polymerization processes is usually more complex; the solution of the mathematical models that describe them can be quite challenging and has been the focus of enormous research efforts. These models may consist of systems of very large dimension: thousands of differential equations, often stiff, which demand special numerical techniques for their solution. In this paper the numerical techniques that can be used to solve this challenging problem are reviewed and contrasted, including weighted residual methods, direct integration, numerical inversion of transformed equations, and lumping methods. Stochastic techniques are also surveyed.     

A Program package using stiff,sparse integration methods for the automatic solution of mass action kinetics equations

A portable program package, MACKSIM, for mass action chemical kinetics simulation, is discussed. As these kinetics are readily expressed in explicit mathematical terms, such a package contains two major and distinct modules, the numerical analysis and the user interface. For the first, MACKSIM uses the latest proven developments incorporating sparse matrix techniques in the backward difference predictor corrector methods originated by Gear for the integration of stiff ordinary differential equations, and thus requires minimal computing time to solve large systems of equations. For the second, the program provides a flexible interface which permits simple specification and variation of reactions, requires no special character input, and has no limit on the number of reactions or species involved other than that imposed by the size of the computer. The technology of these components is discussed briefly, the use of the package for standard reactions is illustrated, and current applications are mentioned.     

Integral tau methods for stiff stochastic chemical systems

Tau leaping methods enable efficient simulation of discrete stochastic chemical systems. Stiff stochastic systems are particularly challenging since implicit methods, which are good for stiffness, result in noninteger states. The occurrence of negative states is also a common problem in tau leaping. In this paper, we introduce the implicit Minkowski-Weyl tau (IMW-τ) methods. Two updating schemes of the IMW-τ methods are presented: implicit Minkowski-Weyl sequential (IMW-S) and implicit Minkowski-Weyl parallel (IMW-P). The main desirable feature of these methods is that they are designed for stiff stochastic systems with molecular copy numbers ranging from small to large and that they produce integer states without rounding. This is accomplished by the use of a split step where the first part is implicit and computes the mean update while the second part is explicit and generates a random update with the mean computed in the first part. We illustrate the IMW-S and IMW-P methods by some numerical examples, and compare them with existing tau methods. For most cases, the IMW-S and IMW-P methods perform favorably.     

Numerical integration of atomic electron density with double exponential formula for density functional calculation

In the preceding study, we reported an application of the double exponential formula to the radial quadrature grid for numerical integration of the radial electron distribution function. Three-type new radial grids with the double exponential transformation were introduced. The performance of radial grids was compared between the double exponential grids and the grids proposed in earlier studies by applying to the electron-counting integrals of noble gas atoms and diatomic molecules including alkali metals, halogens, and transition metals. It was confirmed that the change in accuracy of the quadrature approximation depending on atomic or molecular species is not significant for the double exponential integration schemes rather than the other integration schemes. In the present study, we further investigate the accuracy of the double exponential formula for the electron-counting integrals of all the atoms from H to Kr in the periodic table to elucidate the stable performance of the double exponential radial grids. The electron densities of the atoms are calculated with the Gauss-type orbital basis functions at the B3LYP level. The quadrature accuracy and convergence behavior of numerical integration are compared among the double exponential formula and the formulas proposed by Treutler et al. and by Mura et al. The results reveal that the double exponential radial grids remarkably improve the convergence rate toward high accuracy compared with the previous radial grids, particularly for heavy elements in the 4th period, without fine tuning of the radial grids for each atom.     

Adapted numerical modelling of the Belousov–Zhabotinsky reaction

Adapted numerical schemes for the integration of differential equations generating periodic wavefronts have reported benefits in terms of accuracy and stability. This work is focused on differential equations modelling chemical phenomena which are characterized by an oscillatory dynamics. The adaptation is carried out through the exponential fitting technique, which is specially suitable to follow the apriori known qualitative behavior of the solution. In particular, we have merged this strategy with the information coming from existing theoretical studies and especially the observation of time series. Numerical tests will be provided to show the effectiveness of this problem-oriented approach.     

An equation-free probabilistic steady-state approximation: dynamic application to the stochastic simulation of biochemical reaction networks

Stochastic chemical kinetics more accurately describes the dynamics of "small" chemical systems, such as biological cells. Many real systems contain dynamical stiffness, which causes the exact stochastic simulation algorithm or other kinetic Monte Carlo methods to spend the majority of their time executing frequently occurring reaction events. Previous methods have successfully applied a type of probabilistic steady-state approximation by deriving an evolution equation, such as the chemical master equation, for the relaxed fast dynamics and using the solution of that equation to determine the slow dynamics. However, because the solution of the chemical master equation is limited to small, carefully selected, or linear reaction networks, an alternate equation-free method would be highly useful. We present a probabilistic steady-state approximation that separates the time scales of an arbitrary reaction network, detects the convergence of a marginal distribution to a quasi-steady-state, directly samples the underlying distribution, and uses those samples to accurately predict the state of the system, including the effects of the slow dynamics, at future times. The numerical method produces an accurate solution of both the fast and slow reaction dynamics while, for stiff systems, reducing the computational time by orders of magnitude. The developed theory makes no approximations on the shape or form of the underlying steady-state distribution and only assumes that it is ergodic. We demonstrate the accuracy and efficiency of the method using multiple interesting examples, including a highly nonlinear protein-protein interaction network. The developed theory may be applied to any type of kinetic Monte Carlo simulation to more efficiently simulate dynamically stiff systems, including existing exact, approximate, or hybrid stochastic simulation techniques.     

Modified two-derivative Runge–Kutta methods for the Schrödinger equation

A family of modified two-derivative Runge–Kutta (MTDRK) methods for the integration of the Schrödinger equation are obtained. Two new three-stage and fifth order TDRK methods are derived. The numerical results in the integration of the radial Schrödinger equation with the Woods–Saxon potential are reported to show the high efficiency of our new methods. The results of the error analysis are illustrated by the resonance problem.     

Exponentially fitted multi-derivative linear methods for the resonant state of the Schrödinger equation

A new family of exponentially fitted P-stable one-step linear methods involving several derivatives for the numerical integration of the Schrödinger equation are obtained. Numerical results are reported to show the efficiency and robustness of the new methods specially adapted to the integration of the radial time-independent Schrödinger equation for large energies. Error analysis is carried out and the asymptotic expressions of the local errors for large energies explain the results of the numerical experiments on the resonance problem.     

A family of numerical multistep methods with three distinct schemes: explicit advanced step-point (EAS) methods and the EAS1 approach

In this work, for the first time in an article, we present in a comprehensive way the explicit advanced step-point (EAS) methods. The EAS methods is a family of methods designed for the numerical solution of non-stiff and mildly stiff initial value problems (IVPs) and comprises three distinct schemes: EAS1, EAS2 and EAS3. A thorough theoretical analysis of the EAS family of predictor–corrector methods is presented in terms of their accuracy and stability characteristics and requirements, as well as the rationale for creating the three distinct schemes mentioned above. In this paper we also examine in detail one of the three schemes, the EAS1 methods. EAS1 are assessed for the very first time, are meticulously studied and their superior regions of absolute stability are presented. Furthermore the computational efficiency of EAS1 is examined and comparative numerical results are presented with the use of a variable step, variable order EAS1 code. The numerical results provide good evidence that EAS1 could be seen as superior to the well established Adams methods for the numerical solution of mildly stiff initial value problems.     

A New Eighth-order A-stable Method for Solving Differential Systems Arising in Chemical Reactions

Implicit Runge–Kutta methods are successful algorithms for the numerical solution of stiff differential equations, as they usually appear in chemical reactions. This article describes the construction of a particular implicit method based on internal stages obtained from certain Chebyshev collocation points. The resulting method has algebraic order 8 and A-stability characteristic. An embedding technique using the Runge–Kutta method and a linear multistep one is provided in order to change the step size. Numerical experiments illustrate the behaviour of the new method, showing that it may reach great accuracy and be competitive with other well-known codes.     

Trigonometrically-fitted multi-derivative linear methods for the resonant state of the Schrödinger equation

A family of trigonometrically-fitted multi-derivative linear methods for the numerical integration of the Schrödinger equation are constructed. Numerical results show the efficiency and robustness of the new methods when applied to the radial time-independent Schrödinger equation for large energies. Error analysis is carried out and the asymptotic expressions of the local errors for large energies explain the numerical results in the case of the resonance problem.     

Efficiency of alchemical free energy simulations. I. A practical comparison of the exponential formula, thermodynamic integration, and Bennett's acceptance ratio method

We investigate the relative efficiency of thermodynamic integration, three variants of the exponential formula, also referred to as thermodynamic perturbation, and Bennett's acceptance ratio method to compute relative and absolute solvation free energy differences. Our primary goal is the development of efficient protocols that are robust in practice. We focus on minimizing the number of unphysical intermediate states (λ-states) required for the computation of accurate and precise free energy differences. Several indicators are presented which help decide when additional λ-states are necessary. In all tests Bennett's acceptance ratio method required the least number of λ-states, closely followed by the "double-wide" variant of the exponential formula. Use of the exponential formula in only strict "forward" or "backward" mode was not found to be competitive. Similarly, the performance of thermodynamic integration in terms of efficiency was rather poor. We show that this is caused by the use of the trapezoidal rule as method of numerical quadrature. A systematic study focusing on the optimization of thermodynamic integration is presented in a companion paper.     

Eighth order methods with minimal phase‐lag for accurate computations for the elastic scattering phase‐shift problem

Two new hybrid eighth algebraic order two‐step methods with phase‐lag of order twelve and fourteen are developed for computing elastic scattering phase shifts of the radial Schrödinger equation. Based on these new methods we obtain a new variable‐step procedure for the numerical integration of the Schrödinger equation. Numerical results obtained for the integration of the phase shift problem for the well known case of the Lennard–Jones potential show that these new methods are better than other finite difference methods.