Implicit time integrators for differential algebraic equations in Python

Systems of ordinary differential equations (ODE) is one of the key techniques for describing dynamical systems. It's a natural way of expressing the evolution of a system whether it's a mechanical, electrical, chemical, biological system, etc. From a more contemporary perspective, the solution of ODEs is proving crucial in modeling epidemics. Simulating dynamical systems on discrete computers requires the discretization of ODEs, which is usually performed by numerical integration methods. The requirements for numerical methods increase with the rapid growth of complexity and stiffness of the simulated systems. Therefore, the development of new computationally efficient and stable numerical methods is of certain interest. This work aims to study a more general modelling technique than ODEs: the Differential Algebraic Equations (DAEs). Some mathematical theoretical concepts will be developed to understand and then solve those types of problems. This will come along with how to use them in practice when modelling mechanical, electrical, chemical and control problems. Then, some classical methods (such as BDF or Runge-Kutta) will be turned into DAEs solvers as well as semi-explicit and semi-implicit modifications of the Adams Bashforth Moulton (ABM). Currently, a reference solver is DASSL, a Fortran code based on BDF, through which Linda Petzold has been awarded the J. H. Wilkinson Prize for Numerical Software in 1991. On the other hand, Hairer has studied how Runge-Kutta methods perform compared to BDF in and gave a Fortran implementation of another reference solver based on an implicit Runge-Kutta method called RadauIIa. In this work, semi-explicit and semi-implicit ABM will be compared to these reference solvers. In order to evaluate their performance, six test problems will be proposed and performance will be assessed based on the number of function evaluations. At this point, a numerical stability analysis will have been considered through a 2D Dahlquist test equation. Finally, the implementation of those solvers within the Scipy environment will be explained as well as the different ways of improvements that may lead to boost their performance.