Implicit Runge-Kutta formulas (RKF) are one-step methods with a multi-stage scheme that can be of high order and still retain A-stability. They do not present either the disadvantages of BDF listed above. The main problem associated to RKFs is that the com-putation of the solution requires the solution of a non-linear system of equations, which is excessively expensive to compute when compared to the costs involved in solving multi-step methods. In this section, the computational costs involved in the solution of multi-multi-steps and one-step methods are analyzed and compared.
The general form of a linear multi-step formula (LMF) was introduced in Section 1.1.2 as rewritten by grouping in a constantC the terms calculated in previous solution points, yielding
φ . This is usually done by a modified form of Newton
itera-Numerical integration of ODEs arising in fluid power systems 47
tion. The iterative Newton scheme computes in its m-th iteration the approximation to the solution yn+1 as y(nm+1+1) = y( )nm+1+ Δy( )nm+1, with Δy( )n+m1 given by the linear system
(
Ι−hβ0J) (
Δyn( )m+1)
= −φ( )
y( )nm+1 , (4.3)being J, the Jacobian matrix of the function f evaluated at yn, and I the N N× -dimensional identity matrix. The iteration scheme is repeated until a convergence criterion is achieved. The cost involved in solving a N-dimensional ODE system y′ =j f yj
(
1, ,… yN)
with j = 1…N, by means of a linear multi-step method is then determined from the Newton iteration scheme (4.3). The cost to solve each iteration m consists on:
• One function evaluation of f.
• Evaluation of the Jacobian J and LU-factorization of the iteration matrix(I−hβ0J). The LU-factorization of the iteration matrix is very costly, O (N 3/3) operations. Fortu-nately, the same iteration matrix is employed for all Newton iterations required in one integration step.
The costs in each integration step can be reduced by using the same iteration matrix for a few number of steps. Some BDF codes use this approach in those cases where the Jacobian varies slowly from step to step.
On the other hand, the formulation of one-step implicit RKF is formed with inter-mediate stages ki. The s-stage implicit RKF has the form
1 1
where aij and bi are free-choice parameters, whose values will determine the solution accu-racy of the formula and its numerical stability. The implementation of fully implicit RKF is, by far, more costly than BDF methods. At each integration step, the following non-linear algebraic system of s N× equations must be solved:
( ) ( )
, 1,...,i k =ki− f yn+ ∑h aijkj i= s
φ , (4.5)
being N the dimension of the system y′ = f y( ), and s the number of stages. The non-linear system is solved again with a modified Newton iteration scheme
48 Numerical integration of ODEs arising in fluid power systems
The cost of a Newton iteration applied to a fully implicit RKF method is:
• s function evaluations of f.
• Solution of the linear system (4.7), whose LU-factorization requires O sN
(
( )3/ 3)
op-erations. Usually, the Jacobian J is evaluated once for all the iterations within one inte-gration step, and therefore, only one LU-factorization of the s N× -dimensional itera-tion matrix I−h(A J× ) is required at every integration step.
Such excessive linear algebra costs can be reduced by using Butcher’s technique [Butcher 1976], that exploits the special structure of the iteration matrix in (4.7). By transforming the matrix A-1 to a Jordan canonical form, the LU-factors can be solved in diagonal blocks and therefore the LU-factorization of the iteration matrix is reduced. However, implicit RKF methods are still far from being competitive (in terms of efficiency) to BDF methods, since the latter only requires O (N 3 /3) operations to solve each iteration of the Newton scheme.
One way to reduce significantly the computational costs of implicit RKF is found in the semi-implicit Runge-Kutta formulas (SIRK) [Alexander 1977], also named diagonally-implicit Runge-Kutta formulas (DIRK). SIRK formulas are a particularization of diagonally-implicit RKF in which matrix A is lower triangular (i.e. aij =0for i< j). As a result, the stages ki
in (4.4) can be solved successively i = 1, 2,…, s with only one N-dimensional non-linear system (I−haiiJ) to be solved at each stage. The computational cost per step in DIRK methods is reduced from O ((sN) 3/3) to O (sN 3/3) operations.
A SIRK method with all the diagonal terms of matrix A equal (aii =γ for i = 1,…,s) is called singly diagonally implicit (SDIRK) method:
Numerical integration of ODEs arising in fluid power systems 49
The advantage of having a lower triangular A matrix with equal diagonal terms is that all stages ki can be successively solved by using the same LU-factorization of the iteration matrix(I−hγJ). The dominant costs of a SDIRK method are therefore reduced to s func-tion evaluafunc-tions per iterafunc-tion and one LU-factorizafunc-tion of an N-dimensional iterafunc-tion matrix per integration step, yielding to O (N 3/3) operations.
In the above presented LMF and RK implicit methods, the Newton scheme is iter-ated several times at each integration step until certain convergence criterion is met. A new class of implicit methods, which have the advantage of not using an iteration scheme, is presented next. Rosenbrock formulas introduce the Jacobian term directly into the integra-tion formula instead. Rosenbrock methods were first introduced in [Rosenbrock 1963] and they can be interpreted as the application of only a single Newton iteration at each stage ki
of a SIRK formula. This yields to the following formulation:
1 1
Modified Rosenbrock methods (also known as ROW-methods, or generalized RKFs) can be seen as a generalization of (4.9), since they introduce linear terms of Jkj to the stages ki
for j=1..i. By doing this, more freedom is obtained when establishing order conditions of accuracy and stability properties [Wolfbrandt 1973, Kaps 1979]. The modified Rosenbrock formula is written as
1 1
50 Numerical integration of ODEs arising in fluid power systems
The fact that Rosenbrock methods do not require the solution of non-linear systems make them, potentially, very efficient for the integration of stiff systems of ODEs. The Rosenbrock method in (4.10) is the formula proposed in this thesis for the integration of ODEs arising in fluid power systems. So far, it seems that Rosenbrock formulas have not being used in fluid power systems, although they have proved to be very effective in other applications and in the solution of test equations for low-moderate accuracy requirements.