A comprehensive comparison of ODE solvers for biochemical problems Karol Postawa * , Jerzy Szczygiel , Marek Kula _ zy  nski Department of Chemistry, Wroclaw University of Science and Technology, Wybrze _ ze Wyspia nskiego 27, 50-370, Wroclaw, Poland article info Article history: Received 5 July 2019 Received in revised form 10 March 2020 Accepted 16 April 2020 Available online 25 April 2020 Keywords: Model Simulation Numerical Performance Bioenergy abstract The article is focused on a deep and detailed study on available Ordinary Differential Equations (ODEs) numerical solvers for biochemical and bioprocesses purposes, which are an important part of the renewable energy sector. A wide selection of algorithms is tested - starting from simple, single-step explicit methods, ending with implicit multi-step techniques. These include MATLAB, Python, Cþþ, and C# implementations. The test configuration is an ODEs based model that simulates a biogas pro- duction reactor. The research shows that most of the tested solvers pass the accuracy-test (the difference didn’t exceed 0,07%), however only selected are efficient. Most of Runge-Kutta based methods are slow and require an enormous number of steps (more than 2.5  10 8 ). Only multi-step implicit methods are long term solutions - they provide great accuracy while dealing well with stiff, non-smooth ODEs sets. The best from tested solutions were two MATLAB solvers - ode23s and ode15s, as well as a python solver - the LSODA. The first needed averagely 84,051s of calculation time, and 96465 steps, while ode15s required just 11,529s, performing over 20-times fewer steps. The LSODA is ranked somewhere between them with 18,806s of calculation time and the total number of 23730 steps for tested ODEs set. © 2020 Elsevier Ltd. All rights reserved. 1. Introduction Nowadays computer science is very often supporting decision- making in engineering projects, including bioprocesses like bio- fuel production [1,2]. It’s a very efficient and cost-effective method for the selection and optimization of technology and final products. In the opposite of the experiments in scale, the cost doesn’t in- crease significantly with the number of tries, what’s making it a great tool for testing various possibilities [3]. It needs to be mentioned that a mathematical model can’t replace experiments and pilot-scale tests in all cases, but it’s still a very powerful method for selecting the direction of further research [4]. Bioprocesses usually can be simulated by solving a complex setup of Ordinary Differential Equations (ODE) or in a more advanced approach Differential Algebraic Equations (DAE) [5,6]. The necessity of the differential method is caused by the dynamic nature of this process. The most common and logical way is to describe processes by a set of first-order ODEs supported by some additional simple algebraic equations. By solving it, it is possible to generate the description of all changes in the reactor in respect of time. To create a practice and useable research toll, it is necessary to determinate clear objectives. To be competitive the tool has to provide a platform to test input-output dependencies in a fast way, with adequate detailed feedback. The second requirement is inseparably connected with the complexity of the model [7]. This aspect, justified by accuracy demands, affect calculation perfor- mance. Thus the careful selection of the ODE solving method is a crucial point in creating useful research tools. In literature there’s moderately selection in the topic of ODE solver choice [8,9], how- ever, a more detailed study targeted on bioprocesses is required. Different programming environments offer a wide selection of solvers, which represent various numerical methods for integrating systems of ODEs. In every single case, the selection of the most appropriate should be preceded by a test, since its effectiveness depends on a specific equation set, which describes the process. If the equations of the system contain parts, which could lead to a quick change in the solution, the problem is stiff, and a solver tar- geted to this is required [10]. However the real problem is to determine which ODEs are a part of the class of stiff equations, as there is no precise definition of stiffness, and the solving of stiff ODEs set with incorrect selected numerical methods, despite small * Corresponding author. E-mail address: karol.postawa@pwr.edu.pl (K. Postawa). Contents lists available at ScienceDirect Renewable Energy journal homepage: www.elsevier.com/locate/renene https://doi.org/10.1016/j.renene.2020.04.089 0960-1481/© 2020 Elsevier Ltd. All rights reserved. Renewable Energy 156 (2020) 624e633