Integration Of Ordinary Differential Equations Research Articles

Characterization of Symmetry Properties of First Integrals for Submaximal Linearizable Third-Order ODEs

The relationship between first integrals of submaximal linearizable third-order ordinary differential equations (ODEs) and their symmetries is investigated. We obtain the classifying relations between the symmetries and the first integral for submaximal cases of linear third-order ODEs. It is known that the maximum Lie algebra of the first integral is achieved for the simplest equation and is four-dimensional. We show that for the other two classes they are not unique. We also obtain counting theorems of the symmetry properties of the first integrals for these classes of linear third-order ODEs. For the 5 symmetry class of linear third-order ODEs, the first integrals can have 0, 1, 2, and 3 symmetries, and for the 4 symmetry class of linear third-order ODEs, they are 0, 1, and 2 symmetries, respectively. In the case of submaximal linear higher-order ODEs, we show that their full Lie algebras can be generated by the subalgebras of certain basic integrals.

Local Linearization—Runge–Kutta methods: A class of A-stable explicit integrators for dynamical systems

A new approach for the construction of high order A-stable explicit integrators for ordinary differential equations (ODEs) is theoretically studied. Basically, the integrators are obtained by splitting, at each time step, the solution of the original equation in two parts: the solution of a linear ordinary differential equation plus the solution of an auxiliary ODE. The first one is solved by a Local Linearization scheme in such a way that A-stability is ensured, while the second one can be approximated by any extant scheme, preferably a high order explicit Runge–Kutta scheme. Results on the convergence and dynamical properties of this new class of schemes are given, as well as some hints for their efficient numerical implementation. An specific scheme of this new class is derived in detail, and its performance is compared with some Matlab codes in the integration of a variety of ODEs representing different types of dynamics.

Divergence of solutions of polynomial finite-difference equations

A theorem is proven for kth-order polynomial finite-difference equations that guarantees the divergence of solutions. A ‘basin of divergence’ is characterized and an order of divergence is provided. The basin of divergence is shown to depend on k independent parameters. An unconventional compactification method is used. Applications include the multi-step method in the numerical integration of ordinary differential equations, quadratic equations and the Henon map.

Stabilized overlapping modular time integration of coupled differential-algebraic equations

The modeling of technical systems often leads to differential equations that are coupled by constraints. This results in a set of coupled differential-algebraic equations (DAEs). For efficiency concerns the full system of DAEs is split in several subsystems. The solution proceeds in macro steps. In each macro step, the subsystems are solved separately by different time integration methods, while necessary data from other subsystems is approximated by extrapolation or interpolation. After each macro step, at discrete predefined or adaptively defined communication points, the data is updated. The approximation of data between these communication points (for the current macro step) leads to an additional error in the time integration and may furthermore cause instability. Different from modular time integration of ordinary differential equations (ODEs), this effect cannot be fixed by reducing the stepsize below a small stability bound.In this paper we discuss a stabilization strategy for the stable modular time integration of coupled DAEs. The application of this strategy is illustrated for a simple benchmark problem.

Accelerating multi-dimensional combustion simulations using GPU and hybrid explicit/implicit ODE integration

Simulating multi-dimensional combustion with detailed kinetics often requires solving a large number of ordinary differential equation (ODE) problems at each global time step. In many cases, the ODE integrations account for the bulk of the total wall-clock time for the simulation. This paper introduces CHEMEQ2-GPU – a new explicit stiff ODE solver (based on the existing CHEMEQ2 solver) that exploits the parallel architecture of the modern graphics processing unit (GPU) to accelerate ODE integration in multi-dimensional combustion simulations. We also demonstrate efficient application of the CPU and GPU as co-processors, for further speedup. We describe a hybrid explicit/implicit ODE solver approach that combines the strengths of both solver types running simultaneously on the GPU and CPU, respectively. A dynamic load balancing scheme was used to assign the kinetics ODE integrations over all grid points to either the CPU-based implicit solver DVODE (which is the more efficient solver for highly stiff grid points) or CHEMEQ2-GPU (more efficient for moderately stiff or non-stiff grid points). We demonstrate CHEMEQ2-GPU and the hybrid approach in 3-D simulations of homogeneous charge compression ignition (HCCI) engines. The test cases applied two different n-heptane reaction mechanisms (a large detailed model and a small skeletal model) and three different mesh sizes. Engine simulations were performed using KIVA-CHEMKIN. CHEMEQ2 was about 2–3 times faster than DVODE, with similar prediction accuracy. The CHEMEQ2-GPU speedup relative to CHEMEQ2 increased linearly with the number of grid points for the range of meshes tested in this work. Assuming ideal linear scaling of simulation time with number of processors, the speed of CHEMEQ2-GPU on the Tesla C2050 GPU was equivalent to CHEMEQ2 running on approximately 13 parallel 2.8GHz CPU processors for the finest mesh; and the hybrid solver approach was equivalent to CHEMEQ2 on ∼15 such CPU processors. In summary, CHEMEQ2-GPU provided the additional computing power of 14 parallel CPU processors (for the finest mesh tested) and the hybrid solver approach demonstrated a method to efficiently apply these additional co-processors with existing CPU cores for combustion simulations. CHEMEQ2-GPU scales favorably with the number of grid points and is available by request to the authors. This work presents opportunities for further development, particularly in CPU/GPU load balancing algorithms.

Symmetry Classification of First Integrals for Scalar Linearizable Second‐Order ODEs

Symmetries of the fundamental first integrals for scalar second‐order ordinary differential equations (ODEs) which are linear or linearizable by point transformations have already been obtained. Firstly we show how one can determine the relationship between the symmetries and the first integrals of linear or linearizable scalar ODEs of order two. Secondly, a complete classification of point symmetries of first integrals of such linear ODEs is studied. As a consequence, we provide a counting theorem for the point symmetries of first integrals of scalar linearizable second‐order ODEs. We show that there exists the 0‐, 1‐, 2‐, or 3‐point symmetry cases. It is shown that the maximal algebra case is unique.

Computational Scheme of Ordinary Differential Equations Integration on the Basis of Differential Taylor Transformation with Automatic Step and Order Selection

Computational Scheme of Ordinary Differential Equations Integration on the Basis of Differential Taylor Transformation with Automatic Step and Order Selection

Consistent time integration for transient three‐dimensional contact using the NTS‐method

AbstractThe present work aims to investigate stable time integrators for large deformation contact problems within the framework of the well known node‐to‐segment (NTS)‐method. For this kind of problem, standard time integrators fail to conserve the total energy of the system. To remedy this drawback, we combine a mixed method with the concept of a discrete gradient applied to the aforementioned NTS‐method. In the context of nonlinear elastodynamics stable integrators for ordinary differential equations have been extensively developed and investigated during the last two decades. For contact problems, energy consistent integrators have been developed for the NTS‐method (see e.g. Ref. [1]). (© 2011 Wiley‐VCH Verlag GmbH & Co. KGaA, Weinheim)

Shapes and gravitational fields of rotating two-layer Maclaurin ellipsoids: Application to planets and satellites

The exact solution for the shape and gravitational field of a rotating two-layer Maclaurin ellipsoid of revolution is compared with predictions of the theory of figures up to third order in the small rotational parameter of the theory of figures. An explicit formula is derived for the external gravitational coefficient J 2 of the exact solution. A new approach to the evaluation of the theory of figures based on numerical integration of ordinary differential equations is presented. The classical Radau–Darwin formula is found not to be valid for the rotational parameter ϵ 2 = Ω 2/(2 πG ρ 2) ⩾ 0.17 since the formula then predicts a surface eccentricity that is smaller than the eccentricity of the core–envelope boundary. Interface eccentricity must be smaller than surface eccentricity. In the formula for ϵ 2, Ω is the angular velocity of the two-layer body, ρ 2 is the density of the outer layer, and G is the gravitational constant. For an envelope density of 3000 kg m −3 the failure of the Radau–Darwin formula corresponds to a rotation period of about 3 h. Application of the exact solution and the theory of figures is made to models of Earth, Mars, Uranus, and Neptune. The two-layer model with constant densities in the layers can provide realistic approximations to terrestrial planets and icy outer planet satellites. The two-layer model needs to be generalized to allow for a continuous envelope (outer layer) radial density profile in order to realistically model a gas or ice giant planet.

Motion4D-library extended

Motion4D-library extended

Redesigning combustion modeling algorithms for the Graphics Processing Unit (GPU): Chemical kinetic rate evaluation and ordinary differential equation integration

Detailed modeling of complex combustion kinetics remains challenging and often intractable, due to prohibitive computational costs incurred when solving the associated large kinetic mechanisms. The Graphics Processing Unit (GPU), originally designed for graphics rendering on computer and gaming systems, has recently emerged as a powerful, cost-effective supplement to the Central Processing Unit (CPU) for dramatically accelerating scientific computations. Complex scientific computations are now being performed on the GPU in several research fields, such as quantum chemistry, molecular dynamics, and atmospheric modeling. Here, we present methods for exploiting the highly parallel structure of GPUs for combustion modeling. This paper outlines simple algorithm revisions that can be applied to the majority of existing combustion modeling algorithms for GPU computations. Significant simulation acceleration and predictive capability enhancements were obtained by using these GPU-enhanced algorithms for reaction rate evaluation and in ODE integration. For the demonstrations, we implemented the rate evaluation revisions in CHEMKIN and the ODE integration revisions in DASAC and DVODE and we tested the performance for simulating constant-volume ignition using SENKIN. The simulations using the revised algorithms are more than an order of magnitude faster than the corresponding CPU-only simulations, even for a low-end (double-precision) graphics card. Additionally, the computational time scales less than quadratically with the number of chemical species in the kinetic mechanism when using the GPU, as compared to the super-quadratic scaling normally seen with CPU-only chemical kinetics computations; and the GPU-based revisions do not involve approximations to the detailed kinetics. An analysis of the growth rates of combustion mechanism sizes versus computational capabilities of CPUs and GPUs further reveals the important role that GPUs are expected to play in the future of combustion modeling. Finally, we briefly outline practical steps for effectively transitioning from CPU-only to GPU-enhanced combustion modeling.

Computing the Action of the Matrix Exponential, with an Application to Exponential Integrators

A new algorithm is developed for computing $e^{tA}B$, where A is an $n\times n$ matrix and B is $n\times n_0$ with $n_0 \ll n$. The algorithm works for any A, its computational cost is dominated by the formation of products of A with $n\times n_0$ matrices, and the only input parameter is a backward error tolerance. The algorithm can return a single matrix $e^{tA}B$ or a sequence $e^{t_kA}B$ on an equally spaced grid of points $t_k$. It uses the scaling part of the scaling and squaring method together with a truncated Taylor series approximation to the exponential. It determines the amount of scaling and the Taylor degree using the recent analysis of Al-Mohy and Higham [SIAM J. Matrix Anal. Appl., 31 (2009), pp. 970–989], which provides sharp truncation error bounds expressed in terms of the quantities $\|A^k\|^{1/k}$ for a few values of k, where the norms are estimated using a matrix norm estimator. Shifting and balancing are used as preprocessing steps to reduce the cost of the algorithm. Numerical experiments show that the algorithm performs in a numerically stable fashion across a wide range of problems, and analysis of rounding errors and of the conditioning of the problem provides theoretical support. Experimental comparisons with MATLAB codes based on Krylov subspace, Chebyshev polynomial, and Laguerre polynomial methods show the new algorithm to be sometimes much superior in terms of computational cost and accuracy. An important application of the algorithm is to exponential integrators for ordinary differential equations. It is shown that the sums of the form $\sum_{k=0}^p \varphi_k(A)u_k$ that arise in exponential integrators, where the $\varphi_k$ are related to the exponential function, can be expressed in terms of a single exponential of a matrix of dimension $n+p$ built by augmenting A with additional rows and columns, and the algorithm of this paper can therefore be employed.

Tracing the Pareto frontier in bi-objective optimization problems by ODE techniques

In this paper we present a deterministic method for tracing the Pareto frontier in non-linear bi-objective optimization problems with equality and inequality constraints. We reformulate the bi-objective optimization problem as a parametric single-objective optimization problem with an additional Normalized Normal Equality Constraint (NNEC) similar to the existing Normal Boundary Intersection (NBI) and the Normalized Normal Constraint method (NNC). By computing the so called Defining Initial Value Problem (DIVP) for segments of the Pareto front and solving a continuation problem with a standard integrator for ordinary differential equations (ODE) we can trace the Pareto front. We call the resulting approach ODE NNEC method and demonstrate numerically that it can yield the entire Pareto frontier to high accuracy. Moreover, due to event detection capabilities available for common ODE integrators, changes in the active constraints can be automatically detected. The features of the current algorithm are illustrated for two case studies whose Matlab® code is available as Electronic Supplementary Material to this article.

Application of orthogonal expansions for approximate integration of ordinary differential equations

An approximate method to solve the Cauchy problem for normal and canonical systems of second-order ordinary differential equations is proposed. The method is based on the representation of a solution and its derivative at each integration step in the form of partial sums of series in shifted Chebyshev polynomials of the first kind. A Markov quadrature formula is used to derive the equations for the approximate values of Chebyshev coefficients in the right-hand sides of systems. Some sufficient convergence conditions are obtained for the iterative method solving these equations. Several error estimates for the approximate Chebyshev coefficients and for the solution are given with respect to the integration step size.

Frobenius integrable decompositions for ninth-order partial differential equations of specific polynomial type

Frobenius integrable decompositions are presented for a kind of ninth-order partial differential equations of specific polynomial type. Two classes of such partial differential equations possessing Frobenius integrable decompositions are connected with two rational Bäcklund transformations of dependent variables. The presented partial differential equations are of constant coefficients, and the corresponding Frobenius integrable ordinary differential equations possess higher-order nonlinearity. The proposed method can be also easily extended to the study of partial differential equations with variable coefficients.

Integration of ordinary differential equation with a small parameter via approximate symmetries: Reduction of approximate symmetry algebra to a canonical form

Method of integration of second-order ordinary differential equations with twodimensional Lie symmetry algebras by reducing basic symmetries to canonical forms is extended to second-order equations with a small parameter for their approximate integration using two essential approximate symmetries. Canonical forms of basic operators of corresponding approximate Lie algebras Lr, r = 2, 3, 4, as well as general forms of invariant differential equations and their solutions are presented. The similar problems are also solved for systems of two first-order ordinary differential equations with two approximate symmetries.

Finite difference modelling of rupture propagation with strong velocity-weakening friction

SUMMARY We incorporate rate- and state-dependent friction in explicit finite difference (FD) simulations of mode II dynamic ruptures in elastic media, using the Mimetic Operators Split-Node (MOSN) method, with adjustable order of spatial accuracy (second-, fourth- or mixed-order accurate), including an option that is fourth-order accurate at the fault discontinuity as well as in the elastic volume. At fault points, the rate and state equations combined with the spatially discretized momentum conservation equations form a coupled system of ordinary differential equations (ODEs) for slip velocity and state variable. As a consequence of the rapid damping of velocity perturbations due to the direct effect, this system exhibits numerical stiffness that is inversely proportional to velocity squared. Approximate solutions to this velocity-state system are achieved by two different implicit schemes: (i) a fourth-order Rosenbrock integration of the full system using multiple substeps and (ii) low order integrations (backward Euler and trapezoidal) of the velocity equation, time-staggered with analytic integration of the state equation under the approximation of constant slip velocity over the time step. In assessing the numerical schemes, we use three test problems: ruptures with frictional resistance controlled by (i) a slip evolution law with strong velocity-weakening behaviour at high slip rates, representing thermal weakening due to flash heating of microscopic asperity contacts, (ii) the classic (low-velocity) slip evolution law and (iii) the classic aging evolution law. A convergence analysis is carried out using reference solutions from a spectral boundary integral equation method (BIEM) (a method restricted to homogeneous media, with nominal spectral accuracy in space and second-order accuracy in time for smooth solutions). Errors are measured by root-mean-square differences of fault-plane time histories (slip, slip rate, traction and state). MOSN shows essentially the same convergence rates as BIEM: second-order convergence for slip and state-variable misfits, with slower (but at least first-order) convergence for slip rates and tractions. For a given grid spacing, fourth-order MOSN is as accurate as BIEM for all variables except slip-rate. MOSN-Rosenbrock nominally has fourth-order temporal accuracy for the fault-plane velocity-state ODE integration (compared to lower-order accuracy for the other two MOSN schemes) and therefore provides an important theoretical benchmark. However, it is sensitive to details of the elastic calculation scheme and occasionally its adaptive substepping performs poorly, leading to large excursions from the reference solution. In contrast, MOSN-trapezoidal is robust and reliable, much easier to implement than MOSN-Rosenbrock, and in all cases achieves precision as good as the latter without recourse to substepping. MOSN-Euler has the same advantages as MOSN-trapezoidal, except that its nominal first-order temporal accuracy ultimately leads to larger errors in slip and state variable compared with the higher-order MOSN schemes at sufficiently small grid spacings and time steps.

An updated version of the Motion4D-library

We present an updated version of the Motion4D-library that can be used for the newly developed GeodesicViewer application. New version program summary Program title: Motion4D-library Catalogue identifier: AEEX_v2_0 Program summary URL: http://cpc.cs.qub.ac.uk/summaries/AEEX_v2_0.html Program obtainable from: CPC Program Library, Queen's University, Belfast, N. Ireland Licensing provisions: Standard CPC licence, http://cpc.cs.qub.ac.uk/licence/licence.html No. of lines in distributed program, including test data, etc.: 153 757 No. of bytes in distributed program, including test data, etc.: 5 178 439 Distribution format: tar.gz Programming language: C++ Computer: All platforms with a C++ compiler Operating system: Linux, Unix, Windows RAM: 31 MBytes Catalogue identifier of previous version: AEEX_v1_0 Journal reference of previous version: Comput. Phys. Comm. 180 (2009) 2355 Classification: 1.5 External routines: Gnu Scientific Library (GSL) ( http://www.gnu.org/software/gsl/) Does the new version supersede the previous version?: Yes Nature of problem: Solve geodesic equation, parallel and Fermi–Walker transport in four-dimensional Lorentzian spacetimes. Solution method: Integration of ordinary differential equations. Reasons for new version: To be applicable for the GeodesicViewer (accepted for publication in Comput. Phys. Comm. (COMPHY) 3941, doi:10.1016/j.cpc.2009.10.010 [program AEFP_v1_0]), there were several minor adjustments to be done. Summary of revisions: 1. Calculation of embedding diagrams are improved. 2. Physical units can be used for some metrics. 3. Tests for the constraint equation within the metric classes are slightly modified. 4. New metrics: AlcubierreWarp, GoedelScaled, GoedelScaledCart, Kasner. Running time: The test runs provided with the distribution require only a few seconds to run.

Nonlocal symmetries and integrable ordinary differential equations: ẍ+3xẋ+x3=0 and its generalizations

The equation x+3xx+x3=0 is well known in many areas of mathematics and physics. It possesses the algebra sl(3,R) of Lie point symmetries, hence is equivalent to the equation for a free particle, and both left and right Painleve series. We investigate two higher-dimensional analogs in terms of their symmetry and singularity structures. We find a drastic reduction in symmetry and a loss of some of the singularity properties. From the nonlocal symmetries we are able to determine the complete symmetry group as being represented by a five-dimensional Abelian algebra.

FROBENIUS INTEGRABLE DECOMPOSITIONS FOR TWO CLASSES OF NONLINEAR EVOLUTION EQUATIONS WITH VARIABLE COEFFICIENTS

Frobenius integrable decompositions are introduced for partial differential equations with variable coefficients. Two classes of partial differential equations with variable coefficients are transformed into Frobenius integrable ordinary differential equations. The resulting solutions are illustrated to describe the solution phenomena shared with the KdV and potential KdV equations, the Boussinesq equation and the Camassa–Holm equation with variable coefficients.