A Comparison Between Differential Equation Solver Suites In MATLAB, R, Julia, Python, C, Mathematica, Maple, and Fortran

Parallel numerical integration is a valuable tool used in many applications requiring high-performance numerical solvers, which is of great interest nowadays due to the increasing difficulty and complexity in differential problems. One of the possible approaches to increase the efficiency of ODE solvers is to parallelize recurrent numerical methods, making them more suitable for execution in hardware with natural parallelism, e.g., field-programmable gate arrays (FPGAs) or graphical processing units (GPUs). Some of the simplest and most popular ODE solvers are explicit Runge–Kutta methods. Despite the high implementability and overall simplicity of the Runge–Kutta schemes, recurrent algorithms remain weakly suitable for execution in parallel computers. In this paper, we propose an approach for parallelizing classical explicit Runge–Kutta methods to construct efficient ODE solvers with pipeline architecture. A novel technique to obtain parallel finite-difference models based on Runge–Kutta integration is described. Three test initial value problems are considered to evaluate the properties of the obtained solvers. It is shown that the truncation error of the parallelized Runge–Kutta method does not significantly change after its known recurrent version. A possible speed up in calculations is estimated using Amdahl’s law and is approximately 2.5–3-times. Block diagrams of fixed-point parallel ODE solvers suitable for hardware implementation on FPGA are given.

In the simulation of dynamical processes in economy, social sciences, biology or chemistry, the analyzed values often represent non-negative quantities like the amount of goods or individuals or the density of a chemical or biological species. Such systems are typically described by positive ordinary differential equations (ODEs) that have a non-negative solution for every non-negative initial value. Besides positivity, these processes often are subject to algebraic constraints that result from conservation laws, limitation of resources, or balance conditions and thus the models are differential-algebraic equations (DAEs). In this work, we present conditions under which both these properties, the positivity as well as the algebraic constraints, are preserved in the numerical simulation by Runge---Kutta or multistep discretization methods. Using a decomposition approach, we separate the dynamic and the algebraic equations of a given linear, positive DAE to give positivity preserving conditions for each part separately. For the dynamic part, we generalize the results for positive ODEs to DAEs using the solution representation via Drazin inverses. For the algebraic part, we use the consistency conditions of the discretization method to derive conditions under which this part of the approximation overestimates the exact solution and thus is non-negative. We analyze these conditions for some common Runge---Kutta and multistep methods and observe that for index-1 systems and stiffly accurate Runge---Kutta methods, positivity is conditionally preserved under similar conditions as for ODEs. For higher index problems, however, none of the common methods is suitable.

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.

In recent times, the use of different types of mean in the derivation of explicit Runge-Kutta methods had been on increase. Researchers have explored explicit Runge-Kutta methods derivation by using different types of mean such as geometric mean, harmonic mean, contra-harmonic mean, heronian mean to name but a few; as against the conventional explicit Runge-Kutta methods which was viewed as arithmetic mean. However, despite efforts to improve the derivation of explicit Runge-Kutta methods with use of other types of mean, none has deemed it fit to extend this notion to implicit Runge-Kutta methods. In this article, we present the use of heronian mean as a basis for the construction of implicit Runge-Kutta method in a way of improving the conventional method which is arithmetic mean based. Numerical results was conducted on ordinary differential equations which was compared with the conventional two-stage fourth order implicit Runge-Kutta (IRK4) method and two-stage third order diagonally implicit Runge-Kutta (DIRK3) method. The results presented confirmed that the new scheme performs better than these numerical methods. A better Qualitative properties using Dalquist test equation were established.

In this paper, we study the global truncation error of the linear multistep methods (LMM) in terms of local truncation error of the corresponding Runge-Kutta schemes. The key idea is the representation of LMM with a corresponding Runge-Kutta method. For this, we need to consider the multiple step of a linear multistep method as a single step in the corresponding Runge-Kutta method. Therefore, the global error estimation of a LMM through the Runge-Kutta method will be provided. In this estimation, we do not take into account the effects of roundoff errors. The numerical illustrations show the accuracy and efficiency of the given estimation.

Explicit Runge-Kutta methods (RKMs) are among the most popular classes of formulas for the approximate numerical integration of nonstiff, initial value problems. However, high-order Runge-Kutta methods require more function evaluations per integration step than, for example, Adams methods used in PECE mode, and so, with RKMs, it is expecially important to avoid rejected steps. Steps are often rejected when certain derivatives of the solutions are very large for part of the region of integration. This corresponds, for example, to regions where the solution has a sharp front or, in the limit, some derivative of the solution is discontinuous. In these circumstances the assumption that the local truncation error is changing slowly is invalid, and so any step-choosing algorithm is likely to produce an unacceptable step. In this paper we derive a family of explicit Runge-Kutta formulas. Each formula is very efficient for problems with smooth solution as well as problems having rapidly varying solutions. Each member of this family consists of a fifty-order formula that contains imbedded formulas of all orders 1 through 4. By computing solutions at several different orders, it is possible to detect sharp fronts or discontinuities before all the function evaluations defining the full Runge-Kutta step have been computed. We can then either accpet a lower order solution or abort the step, depending on which course of action seems appropriate. The efficiency of the new algorithm is demonstrated on the DETEST test set as well as on some difficult test problems with sharp fronts or discontinuities.

In the numerical solution of ODEs by implicit time-stepping methods, a system of (nonlinear) equations has to be solved each step. It is common practice to use fixed-point iterations or, in the stiff case, some modified Newton iteration. The convergence rate of such methods depends on the stepsize. Similarly, a stepsize change may force a refactorization of the iteration matrix in the Newton solver. This paper develops new strategies for handling the iterative solution of nonlinear equations in ODE solvers. These include automatic switching between fixed-point and Newton iterations, investigating the optimal convergence rate with respect to total work per unit step, a strategy for when to reevaluate the Jacobian, a strategy for when to refactorize the iteration matrix, coordination with stepsize control. Examples will be given, showing that the new overall strategy works efficiently. In particular, the new strategy admits a restrained stepsize variation without refactorizations, thus permitting the use of a smoother stepsize sequence. The strategy is of equal importance for Runge--Kutta and multistep methods.

Preface. Introduction. 1. Theory of differential equations: an introduction. 1.1 General solvability theory. 1.2 Stability of the initial value problem. 1.3 Direction fields. Problems. 2. Euler's method. 2.1 Euler's method. 2.2 Error analysis of Euler's method. 2.3 Asymptotic error analysis. 2.3.1 Richardson extrapolation. 2.4 Numerical stability. 2.4.1 Rounding error accumulation. Problems. 3. Systems of differential equations. 3.1 Higher order differential equations. 3.2 Numerical methods for systems. Problems. 4. The backward Euler method and the trapezoidal method. 4.1 The backward Euler method. 4.2 The trapezoidal method. Problems. 5. Taylor and Runge-Kutta methods. 5.1 Taylor methods. 5.2 Runge-Kutta methods. 5.3 Convergence, stability, and asymptotic error. 5.4 Runge-Kutta-Fehlberg methods. 5.5 Matlab codes. 5.6 Implicit Runge-Kutta methods. Problems. 6. Multistep methods. 6.1 Adams-Bashforth methods. 6.2 Adams-Moulton methods. 6.3 Computer codes. Problems. 7. General error analysis for multistep methods. 7.1 Truncation error. 7.2 Convergence. 7.3 A general error analysis. Problems. 8. Stiff differential equations. 8.1 The method of lines for a parabolic equation. 8.2 Backward differentiation formulas. 8.3 Stability regions for multistep methods. 8.4 Additional sources of difficulty. 8.5 Solving the finite difference method. 8.6 Computer codes. Problems. 9. Implicit RK methods for stiff differential equations. 9.1 Families of implicit Runge-Kutta methods. 9.2 Stability of Runge-Kutta methods. 9.3 Order reduction. 9.4 Runge-Kutta methods for stiff equations in practice. Problems. 10. Differential algebraic equations. 10.1 Initial conditions and drift. 10.2 DAEs as stiff differential equations. 10.3 Numerical issues: higher index problems. 10.4 Backward differentiation methods for DAEs. 10.5 Runge-Kutta methods for DAEs. 10.6 Index three problems from mechanics. 10.7 Higher index DAEs. Problems. 11. Two-point boundary value problems. 11.1 A finite difference method. 11.2 Nonlinear two-point boundary value problems. Problems. 12. Volterra integral equations. 12.1 Solvability theory. 12.2 Numerical methods. 12.3 Numerical methods - Theory. Problems. Appendix A. Taylor's theorem. Appendix B. Polynomial interpolation. Bibliography. Index.

Many HPC applications that solve differential equations rely on the Runge–Kutta (RK) family of methods for time integration . Among these methods, the fourth-order accurate RK4 scheme is especially popular. This time integration scheme requires applications to evaluate four intermediate stages to take one time step. Depending on the complexity of the problem being solved, the evaluation of these intermediate stages can be computationally expensive. In this paper we develop explicit fourth-order accurate multistep RK methods. The advantage of such methods is that they re-use data from previous time steps, thus requiring fewer intermediate stage evaluations and potentially speeding up applications. We outline a procedure to obtain and tune the method’s coefficients by adjusting their stability regions in an attempt to maximize the size that a time step can take. We validate and evaluate our new methods in the context of numerical relativity applications using the EinsteinToolkit . We believe, however, that these methods and results should generalize to other applications using explicit RK methods.

High-order linear multistep methods with general monotonicity and boundedness properties

In this paper a new method of numerically solving ordinary differential equations is presented. This method is based on the Gaussian numerical integration of different orders. Using two different orders for numerical integration, an adaptive method is derived. Any other numerical solver for ordinary differential equations can be used alongside this method. For instance, Runge–Kutta and Adams–Bashforth–Moulton methods are used together with this new adaptive method. This method is fast, stable, consistent, and suitable for very high accuracies. The accuracy of this method is always higher than the method used alongside with it. Two applications of this method are presented in the field of satellite geodesy. In the first application, for different time periods and sampling rates (increments of time), it is shown that the orbit determined by the new method is—with respect to the Keplerian motion—at least 6 and 25,000,000 times more accurate than, respectively, Runge–Kutta and Adams–Bashforth–Moulton methods of the same degree and absolute tolerance. In the second application, a real orbit propagation problem is discussed for the GRACE satellites. The orbit is propagated by the new numerical solver, using perturbed satellite motion equations up to degree 280. The results are compared with another independent method, the Unscented Kalman Filter. It is shown that the orbit propagated by the numerical solver is approximately 50 times more accurate than the one propagated by the Unscented Kalman Filter approach.

Strong-stability-preserving, Hermite-Birkhoff time-discretization based on k step methods and 8-stage explicit Runge-Kutta methods of order 5 and 4

An explicit Runge-Kutta (RK) method with an extended stability region (ESR) can be constructed based on an RK method with a stability function corresponding to the Taylor series for the exponential function. Such a method modifies the base method's stability region in order either to include (or exclude) some areas of the complex plane or to elongate the stability region along one of the axes. There exist other techniques of constructing RK-ESR methods. This paper discusses their applicability for dynamic simulation of large power systems. It is shown that these methods decrease the accuracy in some areas of the complex plane. In particular, the accuracy degradation relates to low-frequency and/or poorly damped modes; decreasing the time step size in order to improve the accuracy is less efficient. The RK-ESR methods should be considered highly specialized and, in general, are not applicable for dynamic simulation of large power systems.

Parallel implementation for the two-stage SDIRK methods via diagonalization

In the numerical solution of partial differential equations using a method-of-lines approach, the availability of high order spatial discretization schemes motivates the development of sophisticated high order time integration methods. For multiphysics problems with both stiff and nonstiff terms implicit-explicit (IMEX) time stepping methods attempt to combine the lower cost advantage of explicit schemes with the favorable stability properties of implicit schemes. Existing high order IMEX Runge--Kutta or linear multistep methods, however, suffer from accuracy or stability limitations. This work shows that IMEX general linear methods (GLMs) are competitive alternatives to classic IMEX schemes for large problems arising in practice. High order IMEX-GLMs are constructed in the partitioned GLM framework developed earlier by the authors [J. Sci. Comput., 61 (2014), pp. 119--144]. The stability regions of the new schemes are optimized numerically. The resulting IMEX-GLMs have similar stability properties as IMEX Runge--Kutta methods, but they do not suffer from order reduction and are superior in terms of accuracy and efficiency. The new IMEX-GLMs have considerably better stability properties than the IMEX linear multistep methods. Numerical experiments with two- and three-dimensional test problems illustrate the potential of the new schemes to speed up complex applications.