Integration Of Ordinary Differential Equations Research Articles

Accelerating finite-rate chemical kinetics with coprocessors: Comparing vectorization methods on GPUs, MICs, and CPUs

Accurate and efficient methods for solving stiff ordinary differential equations (ODEs) are a critical component of turbulent combustion simulations with finite-rate chemistry. The ODEs governing the chemical kinetics at each mesh point are decoupled by operator-splitting allowing each to be solved concurrently. An efficient ODE solver must then take into account the available thread and instruction-level parallelism of the underlying hardware, especially on many-core coprocessors, as well as the numerical efficiency. A stiff Rosenbrock and a nonstiff Runge–Kutta ODE solver are both implemented using the single instruction, multiple thread (SIMT) and single instruction, multiple data (SIMD) paradigms within OpenCL. Both methods solve multiple ODEs concurrently within the same instruction stream. The performance of these parallel implementations was measured on three chemical kinetic models of increasing size across several multicore and many-core platforms. Two separate benchmarks were conducted to clearly determine any performance advantage offered by either method. The first benchmark measured the run-time of evaluating the right-hand-side source terms in parallel and the second benchmark integrated a series of constant-pressure, homogeneous reactors using the Rosenbrock and Runge–Kutta solvers. The right-hand-side evaluations with SIMD parallelism on the host multicore Xeon CPU and many-core Xeon Phi co-processor performed approximately three times faster than the baseline multithreaded C++ code. The SIMT parallel model on the host and Phi was 13%–35% slower than the baseline while the SIMT model on the NVIDIA Kepler GPU provided approximately the same performance as the SIMD model on the Phi. The runtimes for both ODE solvers decreased significantly with the SIMD implementations on the host CPU (2.5–2.7×) and Xeon Phi coprocessor (4.7–4.9×) compared to the baseline parallel code. The SIMT implementations on the GPU ran 1.5–1.6 times faster than the baseline multithreaded CPU code; however, this was significantly slower than the SIMD versions on the host CPU or the Xeon Phi. The performance difference between the three platforms was attributed to thread divergence caused by the adaptive step-sizes within the ODE integrators. Analysis showed that the wider vector width of the GPU incurs a higher level of divergence than the narrower Sandy Bridge or Xeon Phi. The significant performance improvement provided by the SIMD parallel strategy motivates further research into more ODE solver methods that are both SIMD-friendly and computationally efficient.

A new subdivision algorithm for the flow propagation using polynomial algebras

The Jet Transport method has emerged as a powerful tool for the numerical integration of ordinary differential equations; it uses polynomial expansions to approximate the flow map associated to the differential equation in the neighbourhood of a reference solution. One of the main drawbacks of the method is that the region of accuracy becomes smaller along the integration. In this paper we introduce a procedure to determine a ball covering the set of given initial conditions that keeps the accuracy of the integration within a selected threshold. The paper gives detailed explanations of the algorithm illustrated with some examples of applicability, as well as a comparison with a previous existing method for the same purpose.

GPU-powered model analysis with PySB/cupSODA.

SummaryA major barrier to the practical utilization of large, complex models of biochemical systems is the lack of open-source computational tools to evaluate model behaviors over high-dimensional parameter spaces. This is due to the high computational expense of performing thousands to millions of model simulations required for statistical analysis. To address this need, we have implemented a user-friendly interface between cupSODA, a GPU-powered kinetic simulator, and PySB, a Python-based modeling and simulation framework. For three example models of varying size, we show that for large numbers of simulations PySB/cupSODA achieves order-of-magnitude speedups relative to a CPU-based ordinary differential equation integrator.Availability and implementationThe PySB/cupSODA interface has been integrated into the PySB modeling framework (version 1.4.0), which can be installed from the Python Package Index (PyPI) using a Python package manager such as pip. cupSODA source code and precompiled binaries (Linux, Mac OS/X, Windows) are available at github.com/aresio/cupSODA (requires an Nvidia GPU; developer.nvidia.com/cuda-gpus). Additional information about PySB is available at pysb.org.Supplementary information Supplementary data are available at Bioinformatics online.

Overcoming the lack of kinetic information in biochemical reactions networks

A main aspect in computational modelling of biological systems is the determination of model structure and model parameters. Due to economical and technical reasons, only part of these details are well characterized, while the rest are unknown. To deal with this difficulty, many reverse engineering and parameter estimation methods have been proposed in the literature, however these methods often need an amount of experimental data not always available. In this paper we propose an alternative approach, which overcomes model indetermination solving an Optimization Problem (OP) with an objective function that, similarly to Flux Balance Analysis, is derived from an empirical biological knowledge and does not require large amounts of data. The system behaviour is described by a set of Ordinary Differential Equations (ODE). Model indetermination is resolved selecting time-varying coefficients that maximize/ minimize the objective function at each ODE integration step. Moreover, to facilitate the modelling phase we provide a graphical formalism, based on Petri Nets, which can be used to derive the corresponding ODEs and OP. Finally, the approach is illustrated on a case study focused on cancer metabolism.

О применении метода М.Н. Лагутинского к интегрированию дифференциальных уравнений 1-го порядка. Часть 2. Интегрирование в квадратурах

The method of M.N. Lagutinski (1871-1915) allows to find rational integrals and Darboux polynomials for given differential ring and thus can be used for integration of ordinary differential equations in symbolic form. A realization of Lagutinski method was made under free opensource mathematics software system Sage and will be presented in this article with application for symbolic integration of 1st order differential equations. The second part is devoted to integration of given differential equation d + d with , Q[, ] in quadratures. According to the theorem of M. Singer the problem of integration in quadratures is equivalent to the finding of integrating factor of the form = exp d + d where , Q[, ]. The function can be found as a root of Darboux polynomial for some auxiliary differentiation of the ring Q[, , ]. By Lagutinski method we can find all Darboux polynomials for given differentiation of polynomial ring if degrees of required polynomials are less than given boundary and thus we can find integration factor of the form stated above. The theory and its realization in Sage are tested on numerous examples from standard for Russia text-book by A. F. Filippov.

О применении метода М.Н. Лагутинского к интегрированию дифференциальных уравнений 1-го порядка. Часть 1. Отыскание алгебраических интегралов

The method of M.N. Lagutinski (1871-1915) allows to find rational integrals and Darboux polynomials for given differential ring and thus can be used for integration of ordinary differential equations in symbolic form. A realization of Lagutinski method was made under free opensource mathematics software system Sage and will be presented in this article with application for symbolic integration of 1st order differential equations. In the first part of the article basic concepts of the Lagutinski method is briefly stated for polynomials rings. Then this method is applied to search of algebraic integrated curves for given ordinary differential equations of the form d + d with , Q[, ]. It is shown how the Lagutinski method allows to look for curves of the given order or to prove that there are not such curves. In particular questions about the optimization of computations and integration in micronomials are considered. The theory and its realization in Sage are tested on numerous examples from standard for Russia text-book by A.F. Filippov. Some recommendations for optimization of the Lagutinski method usage are made in the conclusion of the article.

Lie Group Forward Dynamics of Fixed-Wing Aircraft With Singularity-Free Attitude Reconstruction on SO(3)

This paper proposes an approach to formulation and integration of the governing equations for aircraft flight simulation that is based on a Lie group setting, and leads to a nonsingular coordinate-free numerical integration. Dynamical model of an aircraft is formulated in Lie group state space form and integrated by ordinary-differential-equation (ODE)-on-Lie groups Munthe-Kaas (MK) type of integrator. By following such an approach, it is assured that kinematic singularities, which are unavoidable if a three-angles-based rotation parameterization is applied for the whole 3D rotation domain, do not occur in the proposed noncoordinate formulation form. Moreover, in contrast to the quaternion rotation parameterization that imposes additional algebraic constraint and leads to integration of differential-algebraic equations (DAEs) (with necessary algebraic-equation-violation stabilization step), the proposed formulation leads to a nonredundant ODE integration in minimal form. To this end, this approach combines benefits of both traditional approaches to aircraft simulation (i.e., three angles parameterization and quaternions), while at the same time it avoids related drawbacks of the classical models. Besides solving kinematic singularity problem without introducing DAEs, the proposed formulation also exhibits numerical advantages in terms of better accuracy when longer integration steps are applied during simulation and when aircraft motion pattern comprises steady rotational component of its 3D motion. This is due to the fact that a Lie group setting and applied MK integrator determine vehicle orientation on the basis of integration of local (tangent, nonlinear) kinematical differential equations (KDEs) that model process of 3D rotations (i.e., vehicle attitude reconstruction on nonlinear manifold SO(3)) more accurately than “global” KDEs of the classical formulations (that are linear in differential equations part in the case of standard quaternion models).

Experimental characterization of wingtip vortices in the near field using smoke flow visualizations

In order to predict the axial development of the wingtip vortices strength, an accurate theoretical model is required. Several experimental techniques have been used to that end, e.g. PIV or hot-wire anemometry, but they imply a significant cost and effort. For this reason, we have performed experiments using the smoke-wire technique to visualize smoke streaks in six planes perpendicular to the main stream flow direction. Using this visualization technique, we obtained quantitative information regarding the vortex velocity field by means of Batchelor’s model for two chord-based Reynolds numbers, $$Re_c=3.33\times 10^4$$ and $$10^5$$ . Therefore, this theoretical vortex model has been introduced in the integration of ordinary differential equations which describe the temporal evolution of streak lines as function of two parameters: the swirl number, S, and the virtual axial origin, $$\overline{z_0}$$ . We have applied two different procedures to minimize the distance between experimental and theoretical flow patterns: individual curve fitting at six different control planes in the streamwise direction and the global curve fitting which corresponds to all the control planes simultaneously. Both sets of results have been compared with those provided by del Pino et al. (Phys Fluids 23(013):602, 2011b. doi: 10.1063/1.3537791 ), finding good agreement. Finally, we have observed a weak influence of the Reynolds number on the values S and $$\overline{z_0}$$ at low-to-moderate $$Re_c$$ . This experimental technique is proposed as a low cost alternative to characterize wingtip vortices based on flow visualizations.

Derivation of three-derivative Runge-Kutta methods

We introduce an algorithm for a numerical integration of ordinary differential equations in the form of yź = f(y). We extend the two-derivative Runge-Kutta methods (Chan and Tsai, Numer. Algor. 53, 171---194, 2010) to three-derivative Runge-Kutta methods by including the third derivative yźźź=ăź(y)=fźź(y)(f(y),f(y))+fź(y)fź(y)f(y)$y^{\prime \prime \prime }=\hat {g}(y)=f^{\prime \prime }(y)(f(y), f(y))+f^{\prime }(y)f^{\prime }(y)f(y)$. We present an approach based on the algebraic theory of Butcher (Math. Comp. 26, 79---106, 1972) and the źź$\mathcal {B}-$ series theory of Hairer and Wanner (Computing 13, 1---15 (1974)) combined with the methodology of Chan and Chan (Computing 77(3), 237---252, 2006). In this study, special explicit three-derivative Runge-Kutta methods that possess one evaluation of first derivative, one evaluation of second derivative, and many evaluations of third derivative per step are introduced. Methods with stages up to six and of order up to ten are presented. The numerical calculations have been performed on some standard problems and comparisons made with the accessible methods in the literature.

Singularity-free time integration of rotational quaternions using non-redundant ordinary differential equations

A novel ODE time stepping scheme for solving rotational kinematics in terms of unit quaternions is presented in the paper. This scheme inherently respects the unit-length condition without including it explicitly as a constraint equation, as it is common practice. In the standard algorithms, the unit-length condition is included as an additional equation leading to kinematical equations in the form of a system of differential-algebraic equations (DAEs). On the contrary, the proposed method is based on numerical integration of the kinematic relations in terms of the instantaneous rotation vector that form a system of ordinary differential equations (ODEs) on the Lie algebra $\mathit{so}(3)$ of the rotation group $\mathit{SO}(3)$ . This rotation vector defines an incremental rotation (and thus the associated incremental unit quaternion), and the rotation update is determined by the exponential mapping on the quaternion group. Since the kinematic ODE on $\mathit{so}(3)$ can be solved by using any standard (possibly higher-order) ODE integration scheme, the proposed method yields a non-redundant integration algorithm for the rotational kinematics in terms of unit quaternions, avoiding integration of DAE equations. Besides being ‘more elegant’—in the opinion of the authors—this integration procedure also exhibits numerical advantages in terms of better accuracy when longer integration steps are applied during simulation. As presented in the paper, the numerical integration of three non-linear ODEs in terms of the rotation vector as canonical coordinates achieves a higher accuracy compared to integrating the four (linear in ODE part) standard-quaternion DAE system. In summary, this paper solves the long-standing problem of the necessity of imposing the unit-length constraint equation during integration of quaternions, i.e. the need to deal with DAE’s in the context of such kinematical model, which has been a major drawback of using quaternions, and a numerical scheme is presented that also allows for longer integration steps during kinematic reconstruction of large three-dimensional rotations.

Invariants preserving schemes based on explicit Runge–Kutta methods

Numerical integration of ordinary differential equations with some invariants is considered. For such a purpose, certain projection methods have proved its high accuracy and efficiency. Unfortunately, however, sometimes they can exhibit instability. In this paper, a new, highly efficient projection method is proposed based on explicit Runge–Kutta methods. The key there is to employ the idea of the perturbed collocation method, which gives a unified way to incorporate scheme parameters for projection. Numerical experiments confirm the stability of the proposed method.

A partial Lagrangian method for dynamical systems

We develop a new approach termed as a discount free or partial Lagrangian method for construction of first integrals for dynamical systems of ordinary differential equations (ODEs). It is shown how one can utilize the Legendre transformation in a more general setting to provide the equivalence between a current value Hamiltonian and a partial or discount free Lagrangian when it exists. As a consequence, we develop a discount factor free Lagrangian framework to deduce reductions and closed-form solutions via first integrals for ODEs arising from economics by proving three important propositions. The approach is algorithmic and applies to many state variables of the Lagrangian. In order to show its effectiveness, we apply the method to models, one linear and two nonlinear, with one state variable. We obtain new exact solutions for the last model. The discount free Lagrangian naturally arises in economic growth theory and many other economic models when the control variables can be eliminated at the outset which is not always possible in optimal control theory applications of economics. We explain our method with the help of few widely used economic growth models. We point out the difference between this approach and the more general partial Hamiltonian method proposed earlier for a current value Hamiltonian (Naz et al. in Commun Nonlinear Sci Numer Simul 19:3600–3610, 2014) which is applicable in a general setting involving time, state, costate and control variables.

Об одном методе решения нестационарных задач теплопроводности с несимметричными граничными условиями

С использованием дополнительных граничных условий и дополнительной искомой функции в интегральном методе теплового баланса, получено приближенное аналитическое решение нестационарной задачи теплопроводности для бесконечной пластины при несимметричных граничных условиях первого рода. Решение имеет простой вид тригонометрического полинома с коэффициентами, экспоненциально стабилизирующимися во времени. С увеличением числа членов полинома получаемое решение приближается к точному. Введение зависящей от времени дополнительной искомой функции, задаваемой в одной из граничных точек, позволяет свести решение дифференциального уравнения в частных производных к интегрированию обыкновенного дифференциального уравнения. Дополнительные граничные условия находятся в таком виде, чтобы их выполнение искомым решением было эквивалентно выполнению исходного дифференциального уравнения в граничных точках. Показано, что выполнение уравнения лишь в граничных точках приводит к его выполнению и внутри области, минуя интегрирование по пространственной переменной, заменяемого выполнением искомым решением интеграла теплового баланса (осредненного дифференциального уравнения в частных производных). Отсутствие необходимости интегрирования исходного уравнения по пространственной переменной позволяет применять данный метод к решению нелинейных краевых задач с переменными начальными условиями и физическими свойствами среды и др.

Chemistry acceleration with tabulated dynamic adaptive chemistry in a realistic engine with a primary reference fuel

Detailed kinetic reaction mechanisms are necessary for accurate prediction of combustion characteristics such as ignition and emissions in realistic engines. However, the calculation of chemically reacting flows with detailed chemistry is computational expensive due to the large number of species and reactions involved. In this study, a combined approach of dynamic adaptive chemistry (DAC) and in situ adaptive tabulation (ISAT) for efficient chemistry calculations has been implemented into a three-dimensional flow solver to simulate a homogenous charged compression ignition (HCCI) engine with a primary reference fuel (PRF). In the combined method, ISAT speeds up the chemistry calculation by reducing the number of integrations of ordinary differential equations (ODEs) governing chemical kinetics through tabulating and re-using the ODE solutions. At the meantime, DAC accelerates the necessary ODE integrations via the use of locally valid reduced mechanisms, which are obtained using the direct relation graph (DRG) method. The study shows that ISAT–DAC can achieve a speedup factor of about three with accurate prediction of composition and heat release even for the pressure-varying transient engine simulation with significant composition inhomogeneity resulting from wall heat loss. A detailed analysis reveals that the combined method effectively reduces the computational cost through taking advantage of the respective acceleration characteristics of DAC and ISAT at different combustion stages. In the low temperature combustion stage between about 650K and 1000K, even though the reduction in the mechanism size and consequently the ODE integration of the chemical kinetics by DAC is not significant, the combined method can still reach a speed-up factor of more than 100 due to the fact that the tabulated entries can effectively be reused. For the high temperature region, even though the tabulated entries cannot be reused due to the rapid change of the pressure and large composition inhomogeneity resulting from the active combustion and heat loss, DAC can effectively reduce the size of the needed ODE integrations by freezing a significant number of unimportant species and reactions. The study further quantifies the effect of composition inhomogeneity on the computation efficiency in detailed chemistry calculations in realistic engine simulations. For the case considered, the temperature inhomogeneity due to the wall heat loss obviously increases at around top dead center (TDC), reaching a level of 350K difference in temperature within the combustion chamber. Consequently, the overall computational efficiency in chemistry calculation by the combined method has been reduced by 40% compared to the case without heat loss.

Bal: A library for the brute-force analysis of dynamical systems

This paper describes the functionality and usage of bal, a C/C++ library with a Python front-end for the brute-force analysis of continuous-time dynamical systems described by ordinary differential equations (ODEs). bal provides an easy-to-use wrapper for the efficient numerical integration of ODEs and, by detecting intersections of the trajectory with appropriate Poincaré sections, allows to classify the asymptotic trajectory of a dynamical system for bifurcation analysis. Some examples of application are discussed, concerning two-dimensional bifurcation diagrams, Lyapunov exponents and finite-time Lyapunov exponents, basins of attraction, simulation of switching ODE systems, and integration with AUTO, a software package for continuation analysis. Program summaryProgram title: BAL (Library)Catalogue identifier: AEYY_v1_0Program summary URL:http://cpc.cs.qub.ac.uk/summaries/AEYY_v1_0.htmlProgram obtainable from: CPC Program Library, Queen’s University, Belfast, N. IrelandLicensing provisions: Standard CPC licence, http://cpc.cs.qub.ac.uk/licence/licence.htmlNo. of lines in distributed program, including test data, etc.: 65256No. of bytes in distributed program, including test data, etc.: 438133Distribution format: tar.gzProgramming language: C/C++, Python.Computer: PCs and Workstations.Operating system: Any Unix based system such as Linux and Mac OSX.RAM: Variable, depending on the number of dimensions of the dynamical system and duration of the simulation.Classification: 4.3, 4.12.External routines: CVODE, Boost, HDF5Nature of problem: The numerical analysis of continuous-time nonlinear dynamical systems often requires the computation of a large number of solutions of the system, for varying parameter sets, and the subsequent classification of the steady state solution.Solution method: This library provides a convenient interface to the library CVODE for the integration of sets of differential equations and at the same time implements all the machinery required for integrating and classifying the behaviour of the system for varying parameter sets and initial conditions. A high-level Python interface makes the usage of the library amenable also to users who are not familiar with the C/C++ programming languages.

Application of Chebyshev series to integration of ordinary differential equations with rapidly growing solutions

A solution method for systems of ordinary differential equations is described. This method is based on the approximation of right-hand sides by partial sums of shifted Chebyshev series. The coefficients of the series are determined using Markov quadrature formulas. It is shown that the proposed method is more efficient in comparison to the Runge-Kutta and Adams methods when solving differential equations with rapidly growing solutions.

Qualitative analysis of solutions of nonlinear differential equations with the Hukuhara derivative in the space conv (ℝ2)

We suggest a new approach to the qualitative analysis of dynamic properties of solutions of differential equations with the Hukuhara derivative with the use of geometric properties of convex bodies, the comparison method, and the Chaplygin–Wazewski method for approximate integration of ordinary differential equations. We obtain estimates for the area of solutions for a class of nonlinear differential equations with the Hukuhara derivative.

OpenCL parallel integration of ordinary differential equations: Applications in computational dynamics

In many physical problems the use of numerical simulations presents the only path to obtain insight into the behavior and evolution of the system of interest. GPU, CPU and MIC technologies are frequently employed for simulations on computational dynamics and we present results comparing different schemes for the numerical integration of ordinary differential systems (ODEs) in these architectures. The use of adapted methods with low memory storage (Low storage Runge–Kutta methods) gives good results for low precision studies, whereas the Taylor series method provides a powerful technique for high precision. We show how the computation of several dynamics indicators, such as a fast chaos indicator (FLI) or a phase shift indicator in small neuron networks (Central Pattern Generators), can be efficiently computed on these architectures by means of the numerical ODE methods executed through OpenCL. This high computational time reduction allows real-time simulations or generating video media. Program summaryProgram title: OCL-TIDES, OCL-RKCatalogue identifier: AEVW_v1_0Program summary URL:http://cpc.cs.qub.ac.uk/summaries/AEVW_v1_0.htmlProgram obtainable from: CPC Program Library, Queen’s University, Belfast, N. IrelandLicensing provisions: Standard CPC licence, http://cpc.cs.qub.ac.uk/licence/licence.htmlNo. of lines in distributed program, including test data, etc.: 1836No. of bytes in distributed program, including test data, etc.: 10,617Distribution format: tar.gzProgramming language: C and OpenCL.Computer: Any computer with a CPU or a GPU or a Xeon Phi.Operating system: Linux, MacOS, Windows.Has the code been vectorized or parallelized?: Yes, using the parallelization methods from OpenCL.RAM: Problem dependentSupplementary material: A video with a Fast Lyapunov indicator simulation is available (see Appendix A).Classification: 4.3, 6.5.External routines: OpenCL version 1.2Nature of problem: Solution of ODE problems in generic OpenCL environment oriented to large scale independent sets of initial conditions.Solution method: OpenCL parallel integrator based in Runge–Kutta or Taylor Series Method suitable for large sets of independent sets of initial conditions. Both are able to run either on CPU or in a GPU, using the parallelization methods from OpenCL.Running time: Problem dependent, sample tests take less than a minute to run.

Numerical integration of ordinary differential equations with rapidly oscillatory factors

We present a methodology for numerically integrating ordinary differential equations containing rapidly oscillatory terms. This challenge is distinct from that for differential equations which have rapidly oscillatory solutions: here the differential equation itself has the oscillatory terms. Our method generalises Filon quadrature for integrals, and is analogous to integral techniques designed to solve stochastic differential equations and, as such, is applicable to a wide variety of ordinary differential equations with rapidly oscillating factors. The proposed method flexibly achieves varying levels of accuracy depending upon the truncation of the expansion of certain integrals. Users will choose the level of truncation to suit the parameter regime of interest in their numerical integration.

Capturing plasticity effects in overdriven shocks on the finite scale

An ordinary differential equation (ODE) form of the radial return algorithm, which is essentially a Prandtl-Reuss material model, is combined with a strain-rate hardening model to produce an ODE that describes deviatoric stress through a prescribed density rise. An analytical solution is found to the resulting ODE for a specific choice of one of the hardening model’s parameters. That solution is used to prove that if the prescribed density rise is allowed to be infinitely thin, i.e., like a shock in the mathematical sense, the resulting deviatoric stress is still bounded. In other words, the singularity is integrable; integration of the radial return ODE regularizes the infinite strain rate and resulting yield stress in the presence of an ideal shock singularity. The analytical tools developed for this line of thinking are applied to study the variation of deviatoric stress through a nearly shock-like density rise using different density rise profiles, revealing the impact of the shape choice. The tools are also used to compute what rise times are needed to converge upon the correct value of deviatoric stress through a shock; the results indicate that most contemporary hydrocodes cannot be expected to achieve those rise times. A demonstration of connecting the analytical tools to a hydrocode, using surrogate numerical shock shapes, is provided thereby opening the door for using such surrogates to perform sub-grid computations of converged shock behavior for strain-rate hardening materials.