Construction of Rosenbrock–Wanner method Rodas5P and numerical benchmarks within the Julia Differential Equations package

Rosenbrock–Wanner methods for systems of stiff ordinary differential equations are well known since the seventies. They have been continuously developed and are efficient for differential-algebraic equations of index-1, as well. Their disadvantage that the Jacobian matrix has to be updated in every time step becomes more and more obsolete when automatic differentiation is used. Especially the family of Rodas methods has proven to be a standard in the Julia package DifferentialEquations. However, the fifth-order Rodas5 method undergoes order reduction for certain problem classes. Therefore, the goal of this paper is to compute a new set of coefficients for Rodas5 such that this order reduction is reduced. The procedure is similar to the derivation of the methods Rodas4P and Rodas4P2. In addition, it is possible to provide new dense output formulas for Rodas5 and the new method Rodas5P. Numerical tests show that for higher accuracy requirements Rodas5P always belongs to the best methods within the Rodas family.

contains a wide range of solvers for several types of problems. We restrict our considerations to initial value problems of the type M y = f (t, y), y(t 0 ) = y 0 . (1.1) When matrix M is singular, (1.1) is a system of differential-algebraic equations (DAEs), else a system of ordinary differential equations (ODEs). We assume problem (1.1) to be of index not greater than one. For a detailed definition of the index concept see [7]. For solving such problems Rosenbrock-Wanner (ROW) methods are well known, see [4] and [8] for a recent survey. A ROW scheme with stage-number s for problem (1.1) is defined by: 2) h is the stepsize and y 1 is the approximation of the solution y(t 0 +h). The coefficients of the method are γ , α i j , γ i j , and b i define the weights. Moreover, it holds α i = i−1 j=1 α i j and γ i = γ + i−1 j=1 γ i j . ROW methods are alinearly implicit schemes since only a fixed number of s linear systems have to be solved. The index-1 condition guarantees the regularity of the matrix (M − h γ f y ) for sufficiently small stepsizes h > 0, see [4]. A disadvantage compared to implicit Runge-Kutta methods is the requirement of evaluating the Jacobian matrix f y in every timestep.
Within the Julia package DifferentialEquations.jl it is possible to compute the Jacobian by automatic differentiation. Therefore, ROW methods proved to be very efficient for the solution of stiff ODEs and DAEs. For the following analysis we choose ROS3P [10], Rodas3 [21], Rodas4 [4], Rodas4P [24], Rodas4P2 [25] and Rodas5 [3] from the many implemented ROW methods. Moreover, we include Ros3prl2 and Rodas4PR2 [17] which are successors of ROS3P respectively Ros3PL [9] and Rodas4P, but not yet implemented in DifferentialEquations.jl. These schemes are applicable to index-1 DAEs and are stiffly stable. Stiffly stable methods guarantee R(∞) = 0 for the stability function R(z), which is a desired property when solving problem (1.1), see [4,8]. Since in addition all these methods are A-stable they are L-stable as well. The best known method is certainly Rodas4 from Hairer and Wanner [4]. The other schemes considered here were constructed based on this inspiration. The numbers in brackets denote the order of the embedded methods It is well known that ROW methods suffer from order reduction when they are applied to the Prothero-Robinson model, see [15,18,23]: y = λ (y − g(t)) + g (t), y(0) = g(0). (1.4) For a large stiffness parameter |λ| with (λ) << 0 the order may even drop to one. Scholz [23] and Ostermann and Roche [14] derived additional conditions to be fullfilled such that the order is independent on λ. Ostermann and Roche pointed out that the same conditions occur when semi-discretized parabolic partial differential equations (PDEs) are considered. Methods ROS3P, ROS3PL, Ros3prl2, Rodas4P, Rodas4PR2 and Rodas4P2 were developed according to these additional conditions. The letter P stands for "Prothero-Robinson" as well as "parabolic problem". An alternative way to avoid order reduction is considered in [1,2]. Here, the method does not have to fulfill any additional order conditions. In order to achieve a higher stage order, adapted boundary conditions of the partial differential equation are considered in the calculation of the individual stages. The advantage is that every ROW method is suitable for this. The disadvantage is that additional information about the problem to be solved must be included in the stage evaluation. This may become complicated when pipe networks are considered and thus the boundary conditions must take coupling information into account [26]. A numerical comparison of the two approaches is given in Sect. 4.
W methods for ODEs are ROW methods which do not need an exact Jacobian matrix f y in every timestep. Examples for such methods can be found in [17]. Recently, Jax [5] was able to enlarge the class of certain W methods to problems of differential algebraic equations. Unfortunately a huge amount of additional order conditions have to be satiesfied as well. Conditions up to order two are fullfilled by the Rodas4P2 method. Table 1 summarizes the properties of the schemes considered. The order of convergence for some test problems was obtained numerically by the solution with different constant timesteps. The definition of the test problems is given in Sect. 4. Despite the fact that Rodas4 and Rodas5 are very efficient for a couple of typical test problems [4], they show remarkable order reduction for special problems. Rodas5 has some further disadvantages. For simple non-autonomous problems such as y = cos(t), y(0) = 0, errors of this method and its embedded scheme are exactly the same. This leads to a failure of the stepsize control. The embedded method of Rodas5 is not A-stable. In Fig. 1 we can see that the stability domain does not contain the whole left complex half-plane. This may cause stepsize reductions for problems with eigenvalues near the imaginary axis. Moreover, the original literature [3] does not contain a coefficient set for a dense output formula of Rodas5. In the Julia implementation a Hermite interpolation is used which is only applicable to ODE problems.
The aim of this paper therefore is to construct a new coefficient set for Rodas5. It should still have order 5(4) for standard DAE problems of index-1, but its order reduction shown in Table 1 should be restricted to that of Rodas4P2. Moreover, the embedded method should be A-stable and a dense output at least of order p = 4 should be provided. In Sect. 2, all order conditions to be fullfilled by the new method Rodas5P are stated. The construction and the computation of the coefficients of the method is explained in Sect. 3, and finally, in Sect. 4, some numerical benchmarks are given.

Order conditions
The order conditions for Rosenbrock methods applied to index-1 DAEs of type (1.1) were derived by Roche [20]. They are connected to Butcher trees, as shown in Tables 2  and 3. Table 2 lists the conditions up to order p = 5 for ODE problems and Table 3 the additional order conditions up to order p = 5 for index-1 DAE problems, see [3,4].
The following abbreviations are used: The sums in the tables are formed over all possible indices. In Table 4 additional order conditions are defined. Conditions No. 41-44 are given in [11] for problems of type M(y) · y = f (y) with singular matrix M(y). Based on these conditions the method rowdaind2 was derived in [11]. In the special case of index-2 DAEs of type 3) 0 = g(y) (2.4) with a non-singular matrix ( ∂ g ∂ y · ∂ f ∂z ) in the neighborhood of the solution, condition No. 41 guarantees convergence order p = 2. The additonal conditions No. 42-44 lead to order p = 3 for the differential variable y and p = 2 for the algebraic variable z of such index-2 problems.
Conditions No. 45-47 are introduced by Jax [5]. By these conditions at least order p = 2 is achieved for index-1 problems using inexact Jacobian matrices. Method Rodas4P2 has been derived for that purpose, see [25]. When the Jacobian is computed by finite differences, this property might be advantageous. and In order to fulfill conditions No. 48 and 49 in Table 4 all coefficients A i and B i must be zero. As stated in [25], the estimation of the error constant C of the global error in the paper of Scholz [23] is not sharp. It behaves like C = 1 z C 1 for L-stable methods, see [17]. Therefore, for fixed h asymptotically exact results are obtained for |λ| → ∞, but for fixed large stiffness |λ| only order p − 1 is obtained numerically. This can be seen in Table 1. Although the Rodas4P and Rodas4P2 methods satisfy both conditions No. 48 and 49 they only show order p = 3 for the Prothero-Robinson model with large stiffness |λ|, whereas Rodas4PR2 achieves the full order in stiff case.

Construction of Rodas5P
The aim is to construct a method which fullfills all order conditions stated in Tables 2, 3 and 4. Analogously to [3], we choose s = 8 and want to construct a stiffly accurate method with The embedded method with stage numberŝ = 7 is stiffly accurate, too: According to [3,4] we require Therefore, the following 40 coefficients remain to be determined: In the first step we set α 21 = 3γ and β 21 = 0, see [24].  5 , α 52 , α 65 and β 5 as free parameters and try to compute the remaining ones dependent on these. We set β 32 = ( α 3 α 2 ) 2 ( α 3 3 − γ ) and β 3 = 9 2 β 32 and get A 7 = 0,Â 6 = 0, B 6 = 0 from No. 48, 49. 3. We solve the linear system In the next step we get A 5 = 0,Â 4 = 0, B 4 = 0 from the linear system The obtained degree of freedom will be used for fullfilling remaining order conditions in the iteration process later on. 7. Now we can finish the computation of the β-coefficients by solving Values z ∈ C with stabilty function |R(z)| = 1, black for Rodas5P, blue for its embedded method (colour figure online) The stability region of Rodas5P is shown in Fig. 2. It is an A-stable method and due to the stiffly accurate property it is L-stable, too. Next we derive a dense output formula. According to [4] we In order to get a fourth order interpolation conditions No. 1-8 and 18-22 must be fullfilled for the weights b i (τ ). Note that the right hand side of the conditions must be multiplied with τ n , where n is the number of the solid (=black) nodes of the corresponding tree, see [4]. For example, condition No. 21 now reads This condition can be fullfilled by where the first equation for coefficients b i is already true. Thus we have 3 · 13 = 39 linear equations to be satisfied by 3 · s = 24 coefficients. Nevertheless, the solution is possible for the new Rodas5P as well for the known Rodas5 method.
Rodas5P and the new dense output formula for Rodas5 are implemented in the Github repository of the Julia DifferentialEquations package, see https://github.com/SciML/OrdinaryDiffEq.jl. All coefficients of the methods can be found there in particular.

Numerical benchmarks
First we show that the orders given in Table 1 are attained. We solve test problems with known analytical solution y ana (t) by each solver with different numbers of constant stepsizes and compute the numerical errors and orders of convergence. The error is given in the maximum norm at final time: (4.1) The order p is computed by p = log 2 (err 2h /err h ), where err h denotes the error obtained with stepsize h. The following test problems have been treated:     Table  1. Index-1 DAE 1 0 0 0 (2) , t ∈ [2,4] with solution y 1 (t) = ln(t), y 2 (t) = 1 t ln(t). The theoretical orders of convergence are achieved by all methods, see Table 5.

Prothero-Robinson model
with solution y(t) = g(t), see [23,25]. The results are shown in Table 6 and agree with those shown in Table 1. The new method Rodas5P behaves like Rodas4P2 as expected. Computations with different stiffness parameters in the range λ ∈ [10 0 , 10 5 ] show, that only for the method Ros3prl2 the convergence is independent on λ. This includes also mildly stiff problems, where Rodas4PR2 shows order reduction to p = 3.

Parabolic problem
This problem is a slight modification of a similar problem treated in [2]. Function h(x, t) is chosen in order to get the solution u(x, t) = x 3 · e t . The initial values and Dirichlet boundary condition are taken from this solution. Since u(x, t) is cubic in x, x 2 is exact. The numerical results for n x = 1000 space discretization points are given in Table 7. The methods do not achieve the full theoretical order for parabolic problems, shown in Tabel 1. The reason is, that the theory given in [14] assumes linear problems and vanishing boundary conditions. Similar computations with a linear parabolic problem resulted in the full theoretical order. We see further that the embedded method of Rodas5P has nearly order p = 4, too. Nevertheless, the results of the embedded method are slightly worse so that the stepsize control is expected to work. Additionally the results of Rodas5P are compared to the approach chosen in [2]. As proposed there, method GRK4T is applied to problem (4.2). GRK4T [6] is a 4-stage ROW method of order p = 4 for ordinary differential equations. Due to a special choice of its coefficients, it needs only three function evaluations of the righthand side of the ODE system per timestep. Usually, an order reduction to p = 2 would occur when applying it to semi-discretized parabolic problems. This order reduction is prevented by modifications of the boundary conditions in each stage. Technical details can be found in [2]. For comparison, Rodas4 was modified accordingly in addition to GRK4T. Figure 3 shows the results of Rodas5P and the modified methods. For different time stepsizes resulting in different number of function evaluations, the error according to equation (4.1) is plotted. Due to the automatic differentiation for the computation of the Jacobian and the time derivative, only one additional function call was considered respectively. Additonally, the methods were applied with adaptive stepsizes for different tolerances and the error versus elapsed CPU time is shown. While Rodas5P undergoes the small order reduction shown in Table 7, modified GRK4T and Rodas4 have exactly order p = 4. Nevertheless, Rodas5P is more efficient because it has a lower error constant. 4. Index-2 DAE 1 0 0 0 Methods Rodas3, Rodas4, Rodas5 show order reduction to p = 1. All other methods achieve order p = 2, see Table 8. 5. Inexact Jacobian with solution y 1 (t) = sin(t), y 2 (t) = cos(t). Instead of the exact Jacobian we apply J = 0 0 0 2y 2 . According to Jax [5] the derivative of the algebraic equation with respect to the algebraic variable must be exact. We observe the orders shown in Table 1, Rodas5P behaves like Rodas4P2. 6. Dense output We check the dense output formulae of the fourth and fifth order methods via the problem 1 0 0 0 with solution y 1 (t) = y 2 (t) = t n . A method of order p ≥ n should solve this problem exactly within one timestep of size h = 2. After the solution with one timestep we apply the dense output formula to interpolate the solution to times t i = i · h, i = 1, ..., k, h = 2 2 k , k = 1, 2, 3 and compute the resulting maximum error at these timesteps. The numerical errors for different polynomial degrees n of the solution are given in Table 9. Here we can see that the fourth order methods are equipped with dense output formulae of order p = 3, Rodas5 and Rodas5P are able to interpolate with order p = 4.
Next we look at work-precision diagrams and compare the fourth and fifth order methods. In these investigations and in Table 9 Rodas4PR2 was replaced by Rodas4P since no dense output formula is available for Rodas4PR2. The work-precision diagrams are computed for eight different problems by the function WorkPrecisionSet form the Julia package DiffEqDevTools.jl which is part of DifferentialEqua tions.jl. For different tolerances the corresponding computation times and achieved accuracies are evaluated. We show graphs for two different errors: The l 2error is taken from the solution at every timestep and the L 2 -error is taken at 100 evenly  The results for the embedded methods are shown in the lower part of the Table  Table 7 Numerical results (error and order) for problem 3 (parabolic model) Stepsize The results for the embedded methods are shown in the lower part of the Table  spaced points via interpolation. Thus the latter should reflect the error of the dense output formulae. The reference solutions of the problems are computed by Rodas4P2 with tolerances reltol=abstol=10 −14 .

Parabolic problem
We treat again the problem shown in equations (4.2) and Table 7. It turns out that the new method and Rodas4P2 show the best behavior, see Fig. 4. The order reduction of Rodas4 and Rodas5 is clearly visible at the investigated accuracies.

Hyperbolic problem
This problem is discussed in [22,25]. A hyperbolic PDE is discretized by 250 space points. Since the true solution is linear in space variable x, the approximation of the space derivative by first-order finite differences is exact. In Fig. 4 we can see the improved dense output of Rodas5. While the results of Rodas4 and Rodas5 with respect to the l 2 -errors are very similar Rodas5 is much better with respect to the L 2 -errors. In both cases the new method Rodas5P achieves the best numerical results.

Plane pendulum
The pendulum of mass m = 1 and length L can be modeled in cartesian coordinates x(t), y(t) by the equationsẍ = λx, with Lagrange multiplier λ(t) and gravitational constant g. This system is an index-3 DAE which cannot be solved by methods discussed above. By derivation of the algebraic equation with respect to time we achieve the index-2 and index-1 formulation: We

Transistor amplifier
The two-transistor amplifier was intruduced in [19] and further discussed in [12,13]. It consists of eight equations of type (1.1) with index-1. The work-precision diagram shown in Fig. 6 indicates similar behavior for all methods in the l 2 -error. Regarding the L 2 -error Rodas5 and Rodas5P perform best and the improved dense output of Rodas5 is obvious. The new method Rodas5P cannot beat Rodas5 in this case. 5. Water tube problem This example treats the flow of water through 18 tubes which are connected via 13 nodes, see [12,13]. The 49 unknowns of the system are We adapted these equations in order to get a DAE system of index-1. The corresponding results are shown in Fig. 6. It turns out that Rodas5P is slightly more efficient than Rodas5.

Pollution
This is a standard test problem for stiff solvers and contains 20 equations for the chemical reaction part of an air pollution model, see [12,13,27]. The problem is already part of the Julia package SciMLBenchmarks.jl. The results shown in Fig. 7 indicate again that the fifth order methods are preferable.

Photovoltaic network
The new method shall be used in network simulation, see [26]. Therefore we finally simulate a small electric network consisting of a photovoltaic (PV) element, a battery and a consumer with currents i PV (t), i B (t), i C (t). All elements are connected in parallel between two node potentials U 0 (t), U 1 (t). The first node is grounded, at the second node the sum of currents equals zero. The battery is characterized further by its charge q B (t) and an internal voltage u B (t). These seven states are described by equations The third equation describes the consumer that demands a power P(t). It is assumed that P(t) represents a constant power, but it is switched on or off every hour. The  The results for the embedded methods are shown in the lower part of the Table  Fig. 4 Work-precision diagrams for parabolic and hyperbolic problems     Figure 7 shows that in this example the methods Rodas4 and Rodas5P are most suitable.

Conclusion
Based on the construction method for Rodas4P2 a new set of coefficients for Rodas5 could be derived. The new Rodas5P method combines the properties of Rodas5 (high order for standard problems) and Rodas4P or Rodas4P2 (low order reduction for Prothero-Robinson model and parabolic problems). Moreover, it was possible to compute a fourth-order dense output formula for both methods, Rodas5 and Rodas5P.
In all model problems the numerical results of the new method are in the range of the best methods from the class of Rodas schemes studied.
Therefore, Rodas5P can be recommended in the future as a standard method for stiff problems and index-1 DAEs for medium to high accuracy requirements within the Julia package DifferentialEquations.jl. me during implementation of the method in DifferentialEquations.jl and to the two reviewers for valuable comments to improve the paper.
Funding Open Access funding enabled and organized by Projekt DEAL.

Declarations
Conflict of interest The author declare that they have no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.