Parallel implementation for the two-stage SDIRK methods via diagonalization

In the previous chapter, we introduced explicit Runge-Kutta (ERK) methods and demonstrated how they can be implemented as a hierarchy of Python classes. For most ODE systems, replacing the simple forward Euler method with a higher-order ERK method will significantly reduce the number of time steps needed to reach a specified accuracy. Furthermore, it often leads to reduced computation time, since the additional cost per time step is outweighed by the reduced number of steps. However, there exists a class of ODEs known as stiff systems, where all the ERK methods require very small time steps, and any attempt to increase the time step leads to spurious oscillations and possible divergence of the solution. Stiff ODE systems pose a challenge for explicit methods, and they are better addressed by implicit solvers such as implicit Runge-Kutta (IRK) methods. IRK methods are well-suited for stiff problems and can offer substantial reductions in computation time when tackling challenging problems.

In this chapter, the discretization of continuous systems is presented. The explicit and implicit discrete maps are discussed for numerical predictions of continuous systems. Basic discrete schemes are presented which include forward and backward Euler methods, midpoint, and trapezoidal rule method. An introduction to Runge–Kutta methods is presented, and the Taylor series method and second-order Runge–Kutta method are introduced. The explicit Runge–Kutta methods for third and fourth order are systematically presented. The implicit Runge–Kutta methods are discussed based on the polynomial interpolation, which include a generalized implicit Runge–Kutta method, Gauss method, Radau method, and Lotta methods. In addition to one-step methods, implicit and explicit multi-step methods are discussed, including Adams–Bashforth method, Adams–Moulton methods, and explicit and implicit Adams methods.

Abstract—To construct a scheme of implicit Runge-Kutta methods, there are a number of coefficients that must be determined and satisfying consistency properties and Butcher’s simplifying assumptions. In this paper we provide the numerical simulation technique to obtain a scheme of 10th order Implicit Runge-Kutta (IRK10) method. For simulation process, we construct an algorithm to compute all the coefficients involved in the IRK10 scheme. The algorithm is implemented in a language programming (Turbo Pascal) to obtain all the required coefficients in the scheme. To show that our scheme works correctly, we use the scheme to solve Hénon-Heiles system.Keywords—ODEs, 10th order IRK method, numerical technique, Hénon-Heiles system

The dynamics of a differential algebraic equation takes place on a lower dimensional manifold in phase space. Applying a numerical integration scheme, it is natural to ask if and how this geometric property is preserved by the discrete dynamical system. In the index-1 case answers to this question are obtained from the singularly perturbed case treated by Nipp and Stoffer, Numer. Math. 70 (1995), 245–257, for Runge-Kutta methods and in K. Nipp and D. Stoffer, Numer. Math. 74 (1996), 305–323, for linear multistep methods. As main result of this paper it is shown that also for Runge-Kutta methods and linear multistep methods applied to a index-2 problem of Hessenberg form there is a (attractive) invariant manifold for the discrete dynamical system and this manifold is close to the manifold of the differential algebraic equation.

The standard approach to applying IRK methods in the solution of two-point boundary value problems involves the solution of a non-linear system ofn×s equations in order to calculate the stages of the method, wheren is the number of differential equations ands is the number of stages of the implicit Runge-Kutta method. For two-point boundary value problems, we can select a subset of the implicit Runge-Kutta methods that do not require us to solve a non-linear system; the calculation of the stages can be done explicitly, as is the case for explicit Runge-Kutta methods. However, these methods have better stability properties than the explicit Runge-Kutta methods. We have called these new formulas two-point explicit Runge-Kutta (TPERK) methods. Their most important property is that, because their stages can be computed explicity, the solution of a two-point boundary value problem can be computed more efficiently than is possible using an implicit Runge-Kutta method. We have also developed a symmetric subclass of the TPERK methods, called ATPERK methods, which exhibit a number of useful properties.

While implicit Runge–Kutta (RK) methods possess high order accuracy and important stability properties, implementation difficulties and the high expense of solving the coupled algebraic system at each time step are frequently cited as impediments. We present Irksome , a high-level library for manipulating UFL (Unified Form Language) expressions of semidiscrete variational forms to obtain UFL expressions for the coupled Runge–Kutta stage equations at each time step. Irksome works with the Firedrake package to enable the efficient solution of the resulting coupled algebraic systems. Numerical examples confirm the efficacy of the software and our solver techniques for various problems.

Matrix Riccati Differential Equations (MRDEs) \center X ′ = A 21 − X A 11 + A 22 X − X A 12 X , X ( 0 ) = X 0 X’=A_{21}-XA_{11}+A_{22}X-XA_{12}X,\quad X(0)=X_0 , \endcenter where A i j ≡ A i j ( t ) A_{ij}\equiv A_{ij}(t) , appear frequently throughout applied mathematics, science, and engineering. Naturally, the existing conventional Runge-Kutta methods and linear multi-step methods can be adapted to solve MRDEs numerically. Indeed, they have been adapted. There are a few unconventional numerical methods, too, but they are suited more for time-invariant MRDEs than time-varying ones. For stiff MRDEs, existing implicit methods which are preferred to explicit ones require solving nonlinear systems of equations (of possibly much higher dimensions than the original problem itself of, for example, implicit Runge-Kutta methods), and thus they can pose implementation difficulties and also be expensive. In the past, the property of an MRDE which has been most preserved is the symmetry property for a symmetric MRDE; many other crucial properties have been discarded. Besides the symmetry property, our proposed methods also preserve two other important properties — Bilinear Rational Dependence on the initial value, and a Generalized Inverse Relation between an MRDE and its complementary MRDE. By preserving the generalized inverse relation, our methods are accurately able to integrate an MRDE whose solution has singularities. By preserving the property of bilinear dependence on the initial value, our methods also conserve the rank of change to the initial value and a solution’s monotonicity property. Our methods are anadromic ,meaning if an MRDE is integrated by one of our methods from t = τ t\!\!=\!\!\tau to τ + θ \tau \!+\!\theta and then integrated backward from t = τ + θ t\!\!=\!\!\tau \!+\!\theta to τ \tau using the same method, the value at t = τ t\!\!=\!\!\tau is recovered in the absence of rounding errors. This implies that our methods are necessarily of even order of convergence. For time-invariant MRDEs, methods of any even order of convergence are established, while for time-varying MRDEs, methods of order as high as 10 are established; but only methods of order up to 6 are stated in detail. Our methods are semi-implicit, in the sense that there are no nonlinear systems of matrix equations to solve, only linear ones, unlike any pre-existing implicit method. Given the availability of high quality codes for linear matrix equations, our methods can easily be implemented and embedded into any application software package that needs a robust MRDE solver. Numerical examples are presented to support our claims.

Solving delay differential equations using componentwise partitioning by Runge–Kutta method

It is possible to construct fully implicit Runge-Kutta methods like Gaus-Legendre, Radau-IA, Radau-IIA, Lobatto-IIIA, -IIIB, and -IIIC methods of arbitrary high order of convergence. The aim of this paper is to find a new adaptive time stepping for these classes which is based on the embedding technique. Adaptive time step control with embedding is well-known for Runge-Kutta methods, and therefore new embedded methods of order s-1 for the above classes of fully implicit Runge-Kutta methods are constructed. Since these fully implicit methods need the solution of a huge non-linear system of equations different approaches for non-linear equations are discussed and compared. It can be observed that non-linear solvers like the usually used simplified Newton method have a step size restriction if they are applied on higher order methods. We apply our new methods on some lower dimensional ODEs to show that our approach leads to an efficient method.

Variable step-size techniques in continuous Runge-Kutta methods for isospectral dynamical systems

Chapter 8 Numerical analysis of dynamical systems

Implicit Runge–Kutta (IRK) methods are highly effective for solving stiff ordinary differential equations (ODEs) but can be computationally expensive for large-scale problems due to the need of solving coupled algebraic equations at each step. This study improves IRK efficiency by leveraging parallelism to decouple stage computations and reduce communication overhead, specifically we stably decouple a perturbed version of the stage system of equations and recover the exact solution by solving a Sylvester matrix equation with an explicitly known low-rank right-hand side. Two IRK families—symmetric methods and collocation methods—are analyzed, with extensions to nonlinear problems using a simplified Newton method. Implementation details, shared memory parallel code, and numerical examples, particularly for ODEs from spatially discretized PDEs, demonstrate the efficiency of the proposed IRK technique.

This paper proposes applying the two-stage diagonally implicit Runge-Kutta (2S-DIRK) method of numerical integration to the calculation of electromagnetic transients (EMTs) in a power system. The accuracy and the numerical stability of 2S-DIRK are almost the same as those of the trapezoidal method, while 2S-DIRK does not produce sustained numerical oscillation due to a sudden change of an inductor current or a capacitor voltage unlike the trapezoidal method. First, this paper reviews the 2S-DIRK integration scheme and derives the 2S-DIRK formulas of inductors and capacitors for both linear and nonlinear cases. Then, analytical comparisons of 2S-DIRK with the trapezoidal, backward Euler, and Gear-Shichman methods are carried out, and numerical examples which verify the analytical comparisons are shown. Finally, 2S-DIRK is compared with critical damping adjustment (CDA) implemented in electromagnetic transients program (EMTP) for some simulation cases.

We present an efficient library-based numerical method for simulating the Hodgkin-Huxley (HH) neuronal networks. The key components in our numerical method involve (i) a pre-computed high resolution data library which contains typical neuronal trajectories (i.e., the time-courses of membrane potential and gating variables) during the interval of an action potential (spike), thus allowing us to avoid resolving the spikes in detail and to use large numerical time steps for evolving the HH neuron equations; (ii) an algorithm of spike-spike corrections within the groups of strongly coupled neurons to account for spike-spike interactions in a single large time step. By using the library method, we can evolve the HH networks using time steps one order of magnitude larger than the typical time steps used for resolving the trajectories without the library, while achieving comparable resolution in statistical quantifications of the network activity, such as average firing rate, interspike interval distribution, power spectra of voltage traces. Moreover, our large time steps using the library method can break the stability requirement of standard methods (such as Runge-Kutta (RK) methods) for the original dynamics. We compare our library-based method with RK methods, and find that our method can capture very well phase-locked, synchronous, and chaotic dynamics of HH neuronal networks. It is important to point out that, in essence, our library-based HH neuron solver can be viewed as a numerical reduction of the HH neuron to an integrate-and-fire (I&F) neuronal representation that does not sacrifice the gating dynamics (as normally done in the analytical reduction to an I&F neuron).

A myriad of methods, including many optimization and estimation algorithms, require the sensitivities of a dynamic system governed by a set of ordinary differential equations (ODEs). In this paper, the decoupled direct method (DDM) is derived for the calculation of first- and second-order state transition matrices (STMs) using fixed-step-size Runge-Kutta (RK) ODE solvers. Implications for both explicit and implicit RK methods are discussed. Emphasis is placed on the implicit class because of recent interest in the astrodynamics community and because the decoupling of the state propagation and STM calculation eliminates the need to evaluate the Jacobian and Hessian iteratively at each propagation step. Thus, significant computational savings are achievable compared to the traditional propagation of the variational equations, particularly for systems with expensive dynamics. In the current work, two implementations of the DDM are extended both to second-order STMs and to the double-integrator form of the RK equations. One of the strategies presented in this paper, newly derived to include second-order sensitivities, takes advantage of the structure of the system, allowing for the STMs of nonlinear systems to be calculated using efficient linear-algebra techniques rather than iterative methods. This linear-algebra-based approach is shown to be more efficient than the second alternative, and is also more efficient than the propagation of the variational equations for computationally intensive dynamical systems. The effects of important implicit RK customizations, such as variable-fidelity dynamics models and parallelization, are also examined in two examples: a simple, proof-of-concept forced harmonic oscillator, and an Earth-orbiting space object subject to either a low- or high-computational-cost geopotential model.