Runge-Kutta Solution of Initial Value Problems: Methods, Algorithms and Implementation

Abstract: The last decades have seen a strongly increasing use of computers for modeling larger and more complex systems. This has been made possible by a combination of increasing computer power and easy-to-use graphical modeling environments and libraries of components. An important class of systems are described by initial value problems for dynamical systems, and the increasing size and complexity imply that problems are more often either Differential Algebraic Equations (DAE), or stiff Ordinary Differential Equations (ODE). Time-stepping methods for solving DAEs and stiff ODEs have been extensively studied, and the methods have also been implemented in a number of successful codes. Turning a method into a practical implementation requires a number of supporting algorithms, but in spite of their importance we find these to be less well-studied than the time-stepping methods. The supporting algorithms are the key to adaptivity: how to estimate errors, which step-sizes to use, how to start iterations, when to terminate them, when to re-factor Jacobians, etc. The goal of this dissertation is to partially bridge the gap between methods and implementations by studying and improving the algorithms supporting Implicit Runge-Kutta (IRK) methods applied to DAEs and stiff ODEs. Studies of supporting algorithms require both the construction and analysis of simplified models and experimental validation, where we only vary the algorithmic aspect under consideration. The experimental validation has been carried out in the Generic ODE Solving System (Godess), which has been developed in parallel to this thesis and facilitates the necessary means for a strict validation. The dissertation focuses on the iterative process used to compute the intermediate stages in the steps of the IRK methods. The difficulties in the iterative process and the results presented in this dissertation are: The iterative process need accurate predictors to minimize the number of iterations. We construct predictors of orders exceeding the stage order of the method. We also prove that this is only possible by adapting the predictors to the individual stages. The iterations must be terminated after a finite number of iterations. The effects of this are studied in this dissertation and the results include: For algebraically accurate IRK methods starting with an explicit stage we introduce stage derivative reuse (SDR) methods. We prove stiffness independent bounds for the error propagation in the SDR methods. All algebraically accurate IRK methods can employ SDR methods as accurate and inexpensive error estimators. The iterations are in general terminated when the iteration error is below a pre-defined tolerance. We relate the error at the end-point of the step to the errors in the stages. We present a formula relating the tolerance for the iterations to the tolerance for the time-stepping method. The factorization of the Jacobian is often kept for several consecutive steps. An algorithm giving nearly optimal balance between the time spent in the iterations for the stages and the re-factorization and re-computation of the Jacobian is given. The algorithm uses a simple criterion based on the rate of convergence of the iterations. An analysis verifies that the algorithm also has the right asymptotic behavior.