- Email: [email protected]
Multiple shooting using a dichotomically stable integrator for solving differential–algebraic equations
Applied Numerical Mathematics 42 (2002) 117–131 www.elsevier.com/locate/apnum
Multiple shooting using a dichotomically stable integrator for solving differential–algebraic equations Roland England a,∗ , René Lamour b , Jesús López-Estrada c a The Open University, Milton Keynes, UK b Humboldt-University, Berlin, Germany c Universidad Nacional Autónoma de México, México D.F., Mexico
Abstract In previous work by the first author, it has been established that a dichotomically stable discretization is needed when solving a stiff boundary-value problem in ordinary differential equations (ODEs), when sharp boundary layers may occur at each end of the interval. A dichotomically stable implicit Runge–Kutta method, using the 3-stage, fourth-order, Lobatto IIIA formulae, has been implemented in a variable step-size initial-value integrator, which could be used in a multiple-shooting approach. In the case of index-one differential–algebraic equations (DAEs) the use of the Lobatto IIIA formulae has an advantage, over a comparable Gaussian method, that the order is the same for both differential and algebraic variables, and there is no need to treat them separately. The ODE integrator (SYMIRK [R. England, R.M.M. Mattheij, in: Lecture Notes in Math., Vol. 1230, Springer, 1986, pp. 221–234]) has been adapted for the solution of index-one DAEs, and the resulting integrator (SYMDAE) has been inserted into the multiple-shooting code (MSHDAE) previously developed by R. Lamour for differential– algebraic boundary-value problems. The standard version of MSHDAE uses a BDF integrator, which is not dichotomically stable, and for some stiff test problems this fails to integrate across the interval of interest, while the dichotomically stable integrator SYMDAE encounters no difficulty. Indeed, for such problems, the modified version of MSHDAE produces an accurate solution, and within limits imposed by computer word length, the efficiency of the solution process improves with increasing stiffness. For some nonstiff problems, the solution is also entirely satisfactory. 2001 IMACS. Published by Elsevier Science B.V. All rights reserved. Keywords: Differential–algebraic equations; Boundary-value problems; Dichotomic stability; Multiple shooting
1. Introduction Over the last three decades, numerical methods for solving initial-value problems in DAEs have been receiving considerable and increasing attention internationally (see [2,12,14]), the general purpose * Corresponding author.
E-mail address: [email protected] (R. England). 0168-9274/01/$22.00 2001 IMACS. Published by Elsevier Science B.V. All rights reserved. PII: S 0 1 6 8 - 9 2 7 4 ( 0 1 ) 0 0 1 4 5 - 3
118
R. England et al. / Applied Numerical Mathematics 42 (2002) 117–131
integration code DASSL being one of the main achievements for the solution of such problems involving fully implicit index-one DAEs. In contrast, the first results for boundary-value DAE problems have appeared only during the last decade [1,3,15–17], even though these problems arise in important applications such as the description of multibody systems, semiconductor devices, control theory, detonation modelling and chemical kinetics, among others. Of particular importance is the development of a general-purpose code for solving boundary-value problems for fully implicit index-one DAEs, as well as for ODEs. This objective appears already to have been attained using collocation methods [1], although, in the case of a fully implicit system, it is necessary to convert it into a semi-explicit index-two system, thus introducing unnecessary instabilities which have been observed in some numerical results [1,19]. However, such codes are still in a development stage so far as multiple-shooting techniques are concerned [16,17], especially in the case of extremely wellconditioned (stiff) problems [6]. It is important to note that the numerical solution of boundary-value problems in DAEs presents certain difficulties which do not occur for ODEs. It is necessary to formulate the shooting method in such a way that its Jacobian matrix is not singular, as it might be in its most natural formulation; and it is also necessary to determine initial values for each integration which are consistent with the algebraic constraints of the DAE. These difficulties have been overcome by Lamour [16,17], who has developed a multiple-shooting code (MSHDAE) for solving fully implicit index-one differential–algebraic boundaryvalue problems, based on a multiple-shooting scheme which has a nonsingular Jacobian matrix, with dimension comparable to that used in the ODE case. His scheme also determines consistent initial values for each integration. Independently, England and Mattheij [11,22] have established the need for a dichotomically stable discretization when solving a stiff boundary-value problem in ordinary differential equations, with potentially sharp boundary layers at each end of the interval. A dichotomically stable implicit Runge– Kutta method, of order 4, has been implemented in a variable step-size initial-value integrator (SYMIRK), which could be used in a multiple-shooting approach. An explicit third-order method is also required, both as a predictor, and to provide a local error indicator, which has the correct asymptotic behaviour, both for small step-sizes in the boundary layers, and for large step-sizes where the reduced solution is required outside the boundary layers. In the case of differential–algebraic equations, of the form F(x , x, t) = 0 ∈ Rn , it is not possible to reduce the Lobatto IIIA formulae to a single equation, as can be done in the case of ODEs, but instead, a system of 2n equations must be solved for the derivatives Xi+ 1 , Xi+1 , where Xi , Xi are approximations 2
to x(ti ) and x (ti ) respectively. While the Newton iteration matrix for ODEs is of dimension n × n, in the case of DAEs it is larger, of dimension 2n × 2n, but, while the former remains nonsingular, apart from isolated values of h, for ODEs, the latter also remains nonsingular, apart from isolated values of h, for index-one DAEs. The local error indicator would then be obtained in the normal way, by comparing Xi+1 to its predicted value, but in addition to predictors for the variables Xi+1 , Xi+ 1 , it is also necessary to have predictors for 2 the derivatives. The ODE integrator (SYMIRK) has been adapted, for the solution of index-one DAEs, to produce the integrator (SYMDAE). The standard version of MSHDAE uses a BDF integrator, which is not dichotomically stable, and, as is shown in this paper, for some stiff test problems this fails to integrate across the interval of interest, while
R. England et al. / Applied Numerical Mathematics 42 (2002) 117–131
119
the dichotomically stable integrator SYMDAE encounters no difficulty. Indeed, for such problems, the modified version of MSHDAE produces an accurate solution, which does not deteriorate with increasing stiffness within limits imposed by computer word length. For some nonstiff problems, the solution is also entirely satisfactory.
2. Dichotomic stability In [5,8,10] an approach has been developed for solving stiff BVPs, which may often be considered as singularly perturbed problems with some, not necessarily explicit, small parameter ε, which typifies the thickness of possible boundary layers. For example, consider the BVP for n ODEs on the interval [a, b]: dx = x = L(t)x + q(t), dt H0 x(a) + H1 x(b) = c,
(1) (2)
where x(t), q(t), c ∈ Rn and L(t), H0 , H1 ∈ Rn×n . If some rows of the composite matrix [L(t) q(t)] are considerably larger than others, then it may be partitioned in the form: 1 12 1 1 1 11 ˜ q L L (3) [L q] = ε 21 ε 22 ε 2 q˜ L L ˜ with L(t) = O(1), q(t)
= O(1). The solution then consists of a smooth part, together with some fast modes for which x / x = O( 1ε ). In [4] it is shown that, for a well conditioned BVP (1), (2), the fundamental modes of (1) may be expressed in terms of a dichotomic basis, each member of which has either no significant growth or no significant decay. Normally the fast modes remain fast throughout the interval [a, b], and for a well conditioned problem they are negligible outside boundary layers with thickness of order ε. As with any stiff problem, it can be difficult to obtain a solution of the differential equations (1) in the large part of the interval [a, b] where the fast modes are negligible, but still control the stability of the numerical method. The situation is much more difficult than for initial-value problems (IVPs) and in [11] it is shown that it is important to use a discretization which ensures that growing modes of (1) are represented by growing discrete modes, and decaying modes by decaying discrete modes. That is the main feature of a dichotomically stable method. The following definition provides an understanding of the concept for a constant-coefficient system. Definition 1 ([11]). Let the basis solutions of the homogeneous part of an nth-order linear differential equation have components proportional to eλj t , j = 1, . . . , n. Let M h be a discretization method giving an nth-order linear difference equation with corresponding basis solutions of the discrete problem having components with growth {(rjh )i }, where rjh = 1 + hλj + O(h2 ) and i is the step number. Then M h is dichotomically stable on a region R ⊂ C if and only if, for hλj ∈ R, (1) Re(hλj ) 0 ⇒ |rjh | 1, (2) Re(hλj ) 0 ⇒ |rjh | 1.
120
R. England et al. / Applied Numerical Mathematics 42 (2002) 117–131
England and Mattheij [7] investigated methods with this numerical property, and developed [9] a code (SYMIRK) using a fourth-order symmetric implicit (Lobatto IIIA) Runge–Kutta (R–K) formula 1 1 (4) Xi+ 1 = (Xi+1 + Xi ) − h Xi+1 − Xi , 2 2 8 1 2 1 (5) Xi+1 = Xi + h Xi + Xi+ 1 + Xi+1 , 2 6 3 6 where the system of ODEs is written as x = F(x) ∈ Rn . It is important to note that these formulae can be used equally to integrate in either a forward or a backward direction, and, if satisfied exactly, will reproduce the same results in either direction. If this scheme is applied to the model problem (6) x = λ x − f (t) , |λ| 1, then Xi+1 = Xi rλh − λh
1 f 6 i
1 + 23 fi+ 1 + 16 fi+1 − 12 λh(fi+1 − fi ) 2
1 − 12 λh +
1 2 2 λh 12
(7)
where rλh =
1 2 2 1 + 12 λh + 12 λh 1 2 2 1 − 12 λh + 12 λh
is the fourth-order (2, 2)-Padé approximation to the exponential function eλh . Hence ti+1 λh Xi+1 ∼ Xi e − λf (t) dt |λh| → 0 ti
and Xi+1 ∼ Xi − fi + fi+1
|λh| → ∞ ,
(8)
so that the formulae have the correct asymptotic approximation to the solution for |λh| → 0, and, for |λh| → ∞, the numerical solution is almost parallel to the smooth particular solution, which, in the case of the model problem (6), is given to a first approximation by x(t) ∼ f (t) if |λ| 1 in a region where f (t) is smooth. Thus, any errors do not suffer from exponential growth, and the method can work efficiently using large step-sizes, in such a region. The SYMIRK code, outlined in [9], preserves the dichotomic structure of the fundamental solution space by the application of the dichotomically stable formulae (4) and (5). Not only do these formulae make it possible to integrate in either a forward or a backward direction, but they can be used to integrate accurately within the boundary layers, and also to integrate efficiently where the solution behaves smoothly if an appropriate step-size strategy is used.
3. A predictor–corrector pair To develop a variable step-size integrator based on the fourth-order, Lobatto IIIA, Runge–Kutta formulae (4), (5), another method is needed which must have order at least 3, so that the difference
R. England et al. / Applied Numerical Mathematics 42 (2002) 117–131
121
between the two values at the end of each step, which is O(h4 ), may be used as an error indicator. If such a method is to be explicit and is to avoid unstable growth over one step, then it must not use derivative values. Accordingly, an appropriate method is a four-step extrapolation formula, which, in the case of equal step-sizes is p
Xi+1 = 4Xi − 6Xi−1 + 4Xi−2 − Xi−3 .
(9)
If the 3-stage Lobatto IIIA formula is applied to the test problem (6), it produces an approximation over one step which is given by Eq. (7). Hence, for fast modes with large values of |λh|, the numerical solution is given by Eq. (8), and if this numerical approximation is substituted into the four-step extrapolator, assuming for simplicity that the four steps are of the same size, then p
Xi+1 = 4Xi − 6Xi−1 + 4Xi−2 − Xi−3
∼ Xi − fi + (4fi − 6fi−1 + 4fi−2 − fi−3 ) |λh| → ∞ (iv) + O h5 (h → 0) = Xi − fi + fi+1 − h4 fi−1
and so the principal part of the error indicator is essentially due to the interpolation error at the smooth particular solution. The SYMIRK integrator for ODEs uses this extrapolator for both prediction and error estimation.
4. Implementation and step control When solving a nonlinear system of differential equations, written in the form: x = F(x) ∈ Rn it is convenient to eliminate the intermediate value Xi+ 1 from the corrector formulae (4), (5) to produce 2 a single corrector equation:
1 1 2 h 1 (10) Xi+1 − Xi − h Xi+1 + F (Xi+1 + Xi ) − Xi+1 − Xi + Xi = 0, 6 3 2 8 6 )) is the Jacobian matrix of the differential equation, for the unknown vector Xi+1 ∈ Rn . If J(t) = ∂F(x(t ∂x then the Newton iteration matrix for this compressed corrector equation (10) is 1 1 1 (11) I − hJi+1 − hJi+ 1 + h2 Ji+ 1 Ji+1 ∈ Rn×n 2 2 6 3 12 where Ji+1 = J(ti+1 ) and Ji+ 1 = J(ti+ 1 ). 2 2 The exact value of this Newton iteration matrix will rarely be used, and so, in practice, the corrector equation will be solved by some form of modified Newton iteration, which will converge linearly, rather than quadratically as for true Newton iteration. The rate of convergence can be estimated at each iteration, by comparing consecutive norms of the residual of Eq. (10). In the SYMIRK integrator, the most pessimistic, cumulative estimate of the convergence rate is used at each iteration, and on that basis, if it appears that sufficient accuracy will be achieved within a further 3 iterations, then another iteration is taken. If it does not appear that convergence will be achieved within 3 iterations, even with an updated value of the Jacobian matrix J, then the step is abandoned and the step-size is reduced to a level at which convergence may be expected in one iteration. Convergence is considered to be achieved when the norms
122
R. England et al. / Applied Numerical Mathematics 42 (2002) 117–131
of both the residual and the Newton correction are reduced to less than 18 of the tolerance tol, which is specified as a combination of an absolute tolerance, and a tolerance relative to the current value of the solution, component by component. When convergence is achieved, the last Newton correction is not actually applied, so as to avoid the need to reevaluate the function F(x) at ti+1 and ti+ 1 . 2 Having achieved convergence to a satisfactory solution Xi+1 of Eq. (10), it is necessary to evaluate the error indicator which will be used to determine whether or not the step should be accepted, and what step-size should be used for the next step. As indicated earlier, the error indicator Ei+1 will be given as the difference between the final value Xi+1 , which satisfies Eq. (10), and the third-order predicted value p Xi+1 obtained from Eq. (9). On the first three steps of any particular solution, there will be insufficient past data to enable the predictor (9) to be used. For those three steps, it is necessary to use a lower order predictor of the same type, but no useful error indicator is then available until four steps are completed; if the fourth step is rejected, then all four steps must be rejected. Thereafter, the error indicator can be evaluated at every step, and if its norm exceeds the tolerance tol, then the step will be rejected and the step-size reduced by a factor of 0.5. If it is less than 12 tol for five consecutive steps, then the step-size will be increased by a factor
(−1/4) 1 1 if Ei+1 > tol,
Ei+1 / tol 2 32 or 2.0 if it is even smaller. If Ei+1 is between 12 tol and tol, then the step-size will be reduced by a factor
Ei+1 Ei+1
− 0.5 , 1.0 − tol tol which smoothly approximates ( Ei+1 /( 12 tol))(−1/4) for Ei+1 near 12 tol, and 0.5 for Ei+1 near tol. 5. Multiple shooting for DAEs This section will be devoted mainly to the numerical solution of DAEs of the form: F x , x, t = 0 (a t b) subject to boundary conditions G x(a), x(b) = 0,
(12)
(13)
depend only where F(y, x, t), G(xa , xb ) ∈ Rn . Let F and G be sufficiently smooth, the kernel of Fy := ∂F ∂y on t (N(t) := ker(Fy (x , x, t))), and rank Fy =: r = constant. Q(t) denotes a projection onto N(t) (see [12]). Let Q be smooth and P(t) := I − Q(t). By the assumption that the matrix A1 x , x, t := Fy x , x, t + Fx x , x, t Q(t) is nonsingular everywhere we ensure that (12) is a (uniformly) index-one DAE. Contrary to ODEs of dimension n, only r (< n) independent boundary conditions may be fixed. This is guaranteed (see [18]) iff the shooting matrix S := Gxa X(a, a) + Gxb X(b, a)
R. England et al. / Applied Numerical Mathematics 42 (2002) 117–131
satisfies the condition ker S = N(a), im S = M := im(Gxa , Gxb ),
123
(14)
where X(t, s) denotes the fundamental solution of (12). Condition (14) ensures that • (12), (13) is well-posed; • if x denotes a solution of (12), (13), then x is an isolated solution; • considering the correct description of the initial-value conditions in t = s: F x (t), x(t), t = 0, P(s) x(s) − z = 0,
(15)
and denoting the corresponding solution by x(t; s, z), then the IVP (15) has a unique solution at least in a neighbourhood P(s) x (s) − z δ of a solution x . The solution x of an IVP depends on the initial value P(s)z =: u only, i.e., x(t; s, z) = x(t; s, u). Using the subdivision of the interval [a, b]: a = t0 < t1 < · · · < tm−1 < tm = b the boundary and matching conditions, for forward multiple shooting, become G z0 , x(tm ; tm−1 , zm−1 ) = 0, Pi zi − x(ti ; ti−1 , zi−1 ) = 0, i = 1, . . . , m − 1,
(16)
with Pi := P(ti ). For every shooting point ti we need consistent initial values (yi , zi ), which are determined by the solution of Qi yi + Pi vi = 0, F(yi , ui + vi , ti ) = 0, (17) Qi ui = 0, where zi = Pi zi + Qi zi =: ui + vi (see [16,17]). A combination of (16) and (17) yields the shooting operator G(z0 , zm ) + K−1 Q0 u0 , ui+1 − Pi+1 x(ti+1 ; ti , ui ), S(ξ ) := i = 0, . . . , m − 1, Q i yi + Pi vi , F(yi , ui + vi , ti ), with ξ := (u0 , . . . , um−1 , y0 , v0 , . . . , ym−1 , vm−1 ); K must be chosen so that the image of K−1 Q0 has no intersection with Gxa .
124
R. England et al. / Applied Numerical Mathematics 42 (2002) 117–131
To solve the BVP with this shooting method we have to solve the shooting equation S(ξ ) = 0,
(18)
which has a nonsingular Jacobian with respect to ξ , in a neighbourhood of a solution. The essential part of the Jacobian matrix has the same size and also the same block cyclic structure as in the ODE case. The multiple-shooting code is designed in such a way that the directions of integration in every subinterval may be chosen independently. This algorithm is implemented in the code MSHDAE [17], which solves the IVPs (15) by means of backward differentiation formulae (BDFs) in a code modified from that of Skelboe [23]. It may readily be seen that if such an algorithm is applied to a problem such as (1), (2) in ODEs (a special case of DAEs), then the underlying BDF discretization is not dichotomically stable, and may result in an unstable numerical solution if applied to certain types of stiff problem. For DAEs in general, as discussed in [20, 21], a dichotomically stable discretization should also be used, rather than the BDF discretization used in [17], which would be excellent for stiff IVPs but not in the presence of fast growing modes.
6. A dichotomically stable integrator for DAEs The ODE integrator SYMIRK needs some modification in order to solve IVPs for DAEs. The 3-stage Lobatto IIIA method may be written as 5 1 1 X + X 1− X , (19) Xi+ 1 = Xi + h 2 24 i 3 i+ 2 24 i+1 1 2 1 (20) Xi+1 = Xi + h Xi + Xi+ 1 + Xi+1 , 2 6 3 6 and on substitution into the DAE (12) at ti+ 1 and ti+1 , this gives the following system of equations for 2 the unknown derivatives Xi+ 1 , Xi+1 : 2
5 1 1 F Xi+ 1 , Xi + hXi + hXi+ 1 − hXi+1 , ti+ 1 = 0, 2 24 3 24 2 2
1 2 1 F Xi+1 , Xi + hXi + hXi+ 1 + hXi+1 , ti+1 = 0. 6 3 6 2
(21)
It is not possible to eliminate Xi+ 1 , and so reduce this system to a single equation such as (10), and so 2
the larger system in R2n must be solved. The iteration matrix for modified Newton iteration, with respect to the derivative values Xi+ 1 , Xi+1 is 2
Fyi+ 1 + 13 hFxi+ 1 2
2 hFxi+1 3
2
1 − 24 hFxi+ 1 2
Fyi+1 +
1 hFxi+1 6
,
which is nonsingular for an index-one system of DAEs, apart from isolated values of h.
(22)
R. England et al. / Applied Numerical Mathematics 42 (2002) 117–131
125
In addition to predictors for Xi+1 , Xi+ 1 , it is now necessary to have predictors for the derivatives, 2 which, in the case of equal step-sizes, are as follows: p
Xi+1 = 4Xi − 6Xi−1 + 4Xi−2 − Xi−3 , 35 35 21 5 p Xi+ 1 = Xi − Xi−1 + Xi−2 − Xi−3 . 16 16 16 16 2 The heuristics used in the ODE integrator SYMIRK also need some slight modifications for use in the case of DAEs. The use of the Lobatto IIIA formulae (19), (20) has an advantage, over a comparable Gaussian method, that the order of consistency is the same for both the differential and the algebraic variables, and so there is no need to treat them separately. At each iteration, the formulae (5), (4) should be used to determine Xi+1 , Xi+ 1 , and the residuals for (21) should be reduced until they are smaller 2 than the local error tolerance. If the Newton step is (!Xi+ 1 , !Xi+1 ) then, in view of the form of (20), 2
the corresponding change in Xi+1 is given by h( 23 !Xi+ 1 + 16 !Xi+1 ), and so the appropriate norm of 2 this expression should also be reduced until it is smaller than the local error tolerance. The local error indicator would then be obtained in the normal way, by comparing Xi+1 to its predicted value. The ODE integrator SYMIRK has been adapted as outlined above, and the resulting integrator (SYMDAE) has been successfully used to solve a number of index-one DAE problems. The new integrator SYMDAE has also been inserted into the multiple-shooting code (MSHDAE), to facilitate the solution of differential–algebraic BVPs.
7. Test results While a number of tests have been carried out, the problem which has been studied most extensively may be described as follows:
2 t dx 1 1 2 3 − 9x + 12x + 5 x − exp 20 − x 5 − 2 = 0, dt µ
t dx 2 1 2 4 + 12x + 9x + 5 x − 2 exp −20 = 0, dt µ
dx 2 dx 3 t 1 dx 1 (23) 1 − 12 +5 − 100 3x + exp 20 = 0, 9 dt dt dt µ µ
1 2 4 dx dx t 2 dx 2 +9 −5 + 100 3x − exp −20 = 0, 12 dt dt dt µ µ x 3 x 4 − x 5 = 0, subject to boundary conditions: x 1 (0) = 0, where b = µ/20.
x 2 (b) = 0,
x 3 (b) = e,
x 4 (0) = 2,
126
R. England et al. / Applied Numerical Mathematics 42 (2002) 117–131
This problem has a known exact solution given by x 1 (t) x 2 (t) x 3 (t) x 4 (t) x 5 (t)
= = = = =
0, 0, exp(20t/µ), 2 exp(−20t/µ), 2,
(24)
but any perturbation from that solution, at one point, causes a perturbation which grows in both directions away from that point, like exp(±15t). As the stiffness parameter µ is increased, the behaviour of the exact solution is simply stretched out over the interval of increasing length, while the behaviour of the modes exp(±15t) becomes faster and faster in comparison to that of the exact solution involving exp(±20t/µ). Initial-value integrators may have problems during the multiple-shooting process, either because the initial estimates of the solution are not close enough to the exact smooth solution, or because of the accumulation of local truncation errors which are not controlled by the stability of the integrator. In either case, once there is a significant component of the fast growing modes, an integrator is going to need to take very small steps in order to represent it adequately. The idea behind the dichotomically stable integrator (SYMDAE), is that it should be able to control the accumulation of the local truncation errors and the numerical growth of the fast growing modes, so long as the component of those modes remains small, of the same order of magnitude as the local error tolerance for the integration. Results have previously been reported, in [7,9], for a related system of ODEs, using the integrator SYMIRK. Most of the tests on this system of DAEs were carried out using a single shooting interval, with integration in the direction from t = 0 to t = b = µ/20. Initial estimates of the solution were given at t = 0, and to provide a realistic test of the integrators, these were perturbed from the exact solution by an amount of the same order of magnitude as the local error tolerance, typically 10−6 . In the majority of cases, the initial estimates for [x 1 (0), x 2 (0), x 3 (0), x 4 (0)] were [10−4 , 10−4 , 1, 2], though in one case, the first two components were reduced to 10−5 . The initial estimate for x 5 (0) was given as 3, since the initial condition implies that no value of x 5 (0) is in fact used for the integration; however, an initial estimate of some kind is needed for the computation of consistent initial values. For nonstiff problems, with µ in the range from 1 to 10, the two integrators, BDF and SYMDAE, produced very similar numerical results, with final errors of the same order of magnitude (although a little larger than the local error tolerance), even when the initial estimates were perturbed by as much as 10−2 . Usually, a first integration approximated a perturbed solution, while after a single Newton correction to the vector ξ in Eq. (18), a final solution was obtained. For these cases, the dichotomically stable integrator required more function evaluations than the BDF integrator. For slightly more stiff problems, with µ in the range from 100 to 1000, both integrators ran into considerable difficulties, and were generally unable to complete an integration over the shooting interval. To resolve this problem, the interval [0, b] was subdivided into a number m of shooting intervals, as shown in the tables below, thus using a genuine multiple-shooting method. Up to a point, this produced a satisfactory solution, but the condition number of the Jacobian matrix of Eq. (18) was large in these cases, and the shooting algorithm was having increasing difficulties. For larger values of µ, the behaviour of the BDF integrator only deteriorated, as might be expected of a discretization which is not dichotomically stable in the presence of fast growing modes.
R. England et al. / Applied Numerical Mathematics 42 (2002) 117–131
127
However, for µ in the range from 104 to 1018 , excellent results were obtained using SYMDAE, the errors being well within the local error tolerance, with the integrator itself giving increasing accuracy as the stiffness µ increased. Nevertheless, for values of µ greater than approximately 1015 , the integrator SYMDAE began to encounter difficulties. The condition number of the matrix (22) is then of the order 1036 , and, owing to the finite precision of the computer, it becomes impossible to solve Eq. (21) to within the local error tolerance, without a restriction on the maximum step-size used by the integrator. This in turn causes growth of the unwanted fast modes, and makes a satisfactory solution impracticable for µ > 1018 . The principal results are given in Tables 1 and 2. In previous work on ODEs [5,8,10], it has been shown that in addition to a dichotomically stable integrator, it is necessary to develop techniques for localizing a smooth particular solution (a “pathfinder”), and for selecting appropriate shooting intervals. These aspects of the problem have not been dealt with in this paper, and no claim is made to have produced a general purpose code. Nevertheless, multiple shooting has the potential advantage over other methods, such as collocation, that it is not necessary to perform iterations to solve linear problems, such as that given here. The problem described above is a fully implicit index-one problem. In order to solve such a problem using the collocation code COLDAE [1], it is generally necessary to express it as a semi-explicit index-two problem, thereby introducing unnecessary instabilities. A case of such an instability has been Table 1 Results with the BDF integrator ν
µ = 10ν
b = µ/20
m
Iterations
Function calls
Defect
Cond(Jacobian)
0 1 2
1 10 100
0.05 0.50 5
7.7e−9 1.6e−7
2.9 2.6e+3
7720
3.9e−4
1.3e+1
1e3
50
1 3 fails 1 fails
769 1044
3
1 1 1 50 1
Table 2 Results with the dichotomically stable integrator SYMDAE ν
µ = 10ν
b = µ/20
m
Iterations
0 1 2
1 10 1e2
5e−2 0.5 5
1 1 1 10
3 4
1e3 1e4
5 10 15 17 19
1e5 1e10 1e15 1e17 1e19
1 2 fails 11 fails fails 1 1 1 3 3 fails
50 5e2 (x1,2 (0) = 1e−5) 5e3 5e8 5e13 5e15 5e17
1 1 1 1 1 1
Function calls
Defect
Cond
Total
First/last iter.
3286 6041
2032/1216 3544/1216
8.3e−10 4.0e−7
2.9e0 2.6e+3
183916
34654/1976
3.4e−6
4.8e+3
2846 2846 2846 120495 10564909 >479e6
1752/1046 1752/1046 1752/1046 85876/2136 7745524/34860
7.5e−10 2.4e−9 1.2e−8 4.6e−6 2.4e−6 1.6e−8
8.9e0 1.3e0 1.0e0 1.0e0 1.2e0 1.0e0
128
R. England et al. / Applied Numerical Mathematics 42 (2002) 117–131
observed in Ascher and Spiteri [1, Example 3], who suggest that their problem is ill-conditioned. Lamour [19] solves the same problem, using multiple shooting, without any problems of instability. Nevertheless, in order to provide some comparison, this problem has been solved using COLDAE, with a tolerance on the defect of 10−6 , and the results are shown in Table 3. Because of the special structure of this problem, it is also possible to produce an index-one formulation omitting the equation 5 for the unnecessary derivative dxdt , and results are also shown for that reduced formulation. As may be seen, the general purpose code COLDAE produces a solution using fewer function evaluations, and without the restrictions on the stiffness experienced using multiple shooting. On the other hand, as might be expected, multiple shooting, even before developing the necessary techniques for a general purpose code, clearly shows that few or no iterations are required for such a problem, while COLDAE typically requires 5 or 6 iterations. The multiple-shooting code has also been applied to a problem arising in electrical network simulation. It describes a bipolar ring oscillator which was discussed in detail in [13]. Fig. 1 shows the principle network. The current sources J1 and J2 are modelled by one transistor, J3 is constant and an induction is introduced as shown in Fig. 8 in [13]. The problem concerns the autonomous oscillation of the circuit and the results of the integration of an initial value problem up to the oscillation are used as an initial guess. Because of the unknown period of the oscillation, the dimension of the problem is increased from 15 to 16. The initial guess used at the point t = 0 is given by x(0) = [−9.81e−02, −3.24e−01, −6.03e−01, −4.51e−01, −4.49e−01, −8.81e−01, −1.64e−01, −2.41e−02, −4.98e−01, −1.50e+00, −1.01e+00, −4.10e−02, 1.26e−03, 1.26e−03, 3.06e−04, 6.50e+00] and a solution is to be determined in the interval [0, 10−8 ]. The boundary conditions consist of periodic conditions, and a phase condition which fixes the first component of x (see [18]). Shooting over a single interval, both integrators computed the solution, as shown in Table 4. The computed period is 6.50222552e−8. Fig. 2 shows the voltage at node 1 and Fig. 3 the current iss (shown as Jss in Fig. 1). The subroutines used and the complete result file are available from the authors. Table 3 Results with COLDAE ν 0 1 3 5 10 50 100
Index-2 formulation
Index-1 formulation
Final m
Iter.
Function calls
Final m
Iter.
Function calls
20 20 20 20 20 20 20
4 5 6 6 6 6 6
635 745 855 905 905 905 905
20 20 20 20 20 20 20
4 5 5 5 5 5 5
634 744 744 744 644 644 644
Table 4 Results for bipolar ring oscillator problem Integrator
tol
Iterations
Defect
Cond(Jacobian)
BDF SYMDAE
1e−6 1e−6
4 3
3.6e−6 7.5e−7
1.8e+5 1.8e+5
R. England et al. / Applied Numerical Mathematics 42 (2002) 117–131
129
Fig. 1. Bipolar ring oscillator.
Fig. 2. Voltage u1 .
This problem is also fully implicit of index one, and in addition, has nonseparated boundary conditions. In order to put it in a form which could be solved using COLDAE, it would be necessary to express it as a semi-explicit index-two problem with approximately 30 variables. Once again, this would introduce
130
R. England et al. / Applied Numerical Mathematics 42 (2002) 117–131
Fig. 3. Current iss .
instabilities. However, the version of COLDAE which the authors have is unable to solve problems with more than 20 variables, and it is not a straightforward process to increase all the relevant restrictions.
8. Conclusions The numerical solution of boundary-value problems for DAEs presents all the same problems as in the case of ODEs, together with a number of additional difficulties. In particular, using a natural formulation of the multiple-shooting process, the Jacobian matrix of the resulting set of equations is singular. This has been overcome in the code MSHDAE, by combining the computation of consistent initial values, for each integration, with the classical multiple-shooting method. Additionally, the integration method, used for solving each individual initial-value problem, needs to take account of the dichotomic structure of the fundamental modes. A dichotomically stable method for index-one DAEs has been implemented in the integrator SYMDAE, making use of earlier experience in designing such an integrator for ODEs. Numerical results are presented for an example, demonstrating that within a multiple-shooting environment, a dichotomically stable integrator should provide an efficient way of solving very stiff boundary-value DAE problems. A second example shows that it also provides entirely satisfactory results for nonstiff and nonlinear problems. The multiple-shooting code MSHDAE does not yet have all the components which would be required of a general purpose code, some of which would improve its efficiency. Nevertheless, it has been shown here that, with an appropriate integrator, it is capable of obtaining good solutions of fully implicit and nonlinear index-one problems, even with nonseparated boundary conditions. Such problems usually require more iterations using the collocation code COLDAE, and can suffer from an increase both in the number of variables and, more importantly, the index, which can be the cause of instabilities.
R. England et al. / Applied Numerical Mathematics 42 (2002) 117–131
131
References [1] U.M. Ascher, R.J. Spiteri, Collocation software for boundary value differential–algebraic equations, SIAM J. Sci. Comput. 15 (1994) 938–952. [2] K.E. Brenan, S.L. Campbell, L.R. Petzold, Numerical Solution of Initial-Value Problems in Differential–Algebraic Equations, 2nd edn., SIAM, Philadelphia, PA, 1996. [3] K.D. Clark, L.R. Petzold, Numerical solution of boundary value problems in differential–algebraic systems, SIAM J. Sci. Statist. Comput. 10 (5) (1989) 915–936. [4] F.R. de Hoog, R.M.M. Mattheij, On dichotomy and well conditioning in BVP, SIAM J. Numer. Anal. 24 (1) (1987) 89–105. [5] E. Dueñas, R. England, J. López-Estrada, Multiple shooting with dichotomically stable formulae for linear boundary-value problems, Comput. Math. Appl. 38 (9–10) (1999) 143–159. [6] R. England, R. Lamour, Dichotomic stability and multiple shooting for numerical boundary-value problems in index-one differential–algebraic equations, Z. Angew. Math. Mech. 76 (Supp. 1) (1996) 87–90. [7] R. England, R.M.M. Mattheij, Discretizations with dichotomic stability for two-point boundary value problems, in: U.M. Ascher, R.D. Russell (Eds.), Numerical Boundary Value ODEs, Progress in Scientific Computing, Vol. 5, Birkhäuser, Boston, 1985, pp. 91–106. [8] R. England, R.M.M. Mattheij, A stable sequential marching approach for two-point boundary value problems, Presented at the Workshop on Numerical Analysis and its Applications, Cuarto Coloquio de Matemáticas del Centro de Investigación y de Estudios Avanzados del Instituto Politécnico Nacional, Taxco, Mexico, 1985. [9] R. England, R.M.M. Mattheij, Sequential step control for integration of two-point boundary value problems, in: J.P. Hennart (Ed.), Numerical Analysis Proceedings, Guanajuato, Lecture Notes in Math., Vol. 1230, Springer, Berlin, 1986, pp. 221–234. [10] R. England, R.M.M. Mattheij, A sequential algorithm for solving boundary value problems with large Lipschitz constants, Rep. Appl. Numer. Anal. 87-15, Eindhoven University of Technology, Eindhoven, 1987. [11] R. England, R.M.M. Mattheij, Boundary value problems and dichotomic stability, SIAM J. Numer. Anal. 25 (5) (1988) 1037–1054. [12] E. Griepentrog, R. März, Differential–Algebraic Equations and their Numerical Treatment, Teubner, Leipzig, 1986. [13] M. Günter, U. Feldmann, CAD-based electric-circuit modeling in industry: I. Mathematical structure and index of network equations, Surveys Math. Indust. 8 (1999) 97–129. [14] E. Hairer, G. Wanner, Solving Ordinary Differential Equations II: Stiff and Differential–Algebraic Problems, Springer, New York, 1991. [15] M. Kiehl, Sensitivity analysis of ODEs and DAEs—Theory and implementation guide, Optim. Methods Software 10 (6) (1999) 803–821. [16] R. Lamour, Shooting methods for transferable DAEs—Solution of shooting equation, in: Numerical Analysis and Mathematical Modelling, Vol. 24, Banach Center Publications, Warsaw, 1990, pp. 233–240. [17] R. Lamour, A well-posed shooting method for transferable DAEs, Numer. Math. 59 (1991) 815–829. [18] R. Lamour, Oscillation in DAEs, Preprint Nr. 272, Humboldt-University of Berlin, Department of Mathematics, Berlin, 1991. [19] R. Lamour, A shooting method for fully implicit index-2 differential–algebraic equations, SIAM J. Sci. Comput. 18 (1) (1997) 94–114. [20] M. Lentini, R. März, The conditioning of boundary value problems in transferable differential–algebraic equations, SIAM J. Numer. Anal. 27 (1990) 1001–1015. [21] M. Lentini, R. März, Conditioning and dichotomy in differential algebraic equations, SIAM J. Numer. Anal. 27 (1990) 1519–1526. [22] R.M.M. Mattheij, Decoupling and stability of algorithms for boundary value problems, SIAM Rev. 27 (1) (1985) 1–44. [23] S. Skelboe, INTGR for integration of stiff systems of ordinary differential equations, Instituttet for Teleteknik, Report IT 9, Lyngby, 1977.
Multiple shooting using a dichotomically stable integrator for solving differential–algebraic equations Roland England a,∗ , René Lamour b , Jesús López-Estrada c a The Open University, Milton Keynes, UK b Humboldt-University, Berlin, Germany c Universidad Nacional Autónoma de México, México D.F., Mexico
Abstract In previous work by the first author, it has been established that a dichotomically stable discretization is needed when solving a stiff boundary-value problem in ordinary differential equations (ODEs), when sharp boundary layers may occur at each end of the interval. A dichotomically stable implicit Runge–Kutta method, using the 3-stage, fourth-order, Lobatto IIIA formulae, has been implemented in a variable step-size initial-value integrator, which could be used in a multiple-shooting approach. In the case of index-one differential–algebraic equations (DAEs) the use of the Lobatto IIIA formulae has an advantage, over a comparable Gaussian method, that the order is the same for both differential and algebraic variables, and there is no need to treat them separately. The ODE integrator (SYMIRK [R. England, R.M.M. Mattheij, in: Lecture Notes in Math., Vol. 1230, Springer, 1986, pp. 221–234]) has been adapted for the solution of index-one DAEs, and the resulting integrator (SYMDAE) has been inserted into the multiple-shooting code (MSHDAE) previously developed by R. Lamour for differential– algebraic boundary-value problems. The standard version of MSHDAE uses a BDF integrator, which is not dichotomically stable, and for some stiff test problems this fails to integrate across the interval of interest, while the dichotomically stable integrator SYMDAE encounters no difficulty. Indeed, for such problems, the modified version of MSHDAE produces an accurate solution, and within limits imposed by computer word length, the efficiency of the solution process improves with increasing stiffness. For some nonstiff problems, the solution is also entirely satisfactory. 2001 IMACS. Published by Elsevier Science B.V. All rights reserved. Keywords: Differential–algebraic equations; Boundary-value problems; Dichotomic stability; Multiple shooting
1. Introduction Over the last three decades, numerical methods for solving initial-value problems in DAEs have been receiving considerable and increasing attention internationally (see [2,12,14]), the general purpose * Corresponding author.
E-mail address: [email protected] (R. England). 0168-9274/01/$22.00 2001 IMACS. Published by Elsevier Science B.V. All rights reserved. PII: S 0 1 6 8 - 9 2 7 4 ( 0 1 ) 0 0 1 4 5 - 3
118
R. England et al. / Applied Numerical Mathematics 42 (2002) 117–131
integration code DASSL being one of the main achievements for the solution of such problems involving fully implicit index-one DAEs. In contrast, the first results for boundary-value DAE problems have appeared only during the last decade [1,3,15–17], even though these problems arise in important applications such as the description of multibody systems, semiconductor devices, control theory, detonation modelling and chemical kinetics, among others. Of particular importance is the development of a general-purpose code for solving boundary-value problems for fully implicit index-one DAEs, as well as for ODEs. This objective appears already to have been attained using collocation methods [1], although, in the case of a fully implicit system, it is necessary to convert it into a semi-explicit index-two system, thus introducing unnecessary instabilities which have been observed in some numerical results [1,19]. However, such codes are still in a development stage so far as multiple-shooting techniques are concerned [16,17], especially in the case of extremely wellconditioned (stiff) problems [6]. It is important to note that the numerical solution of boundary-value problems in DAEs presents certain difficulties which do not occur for ODEs. It is necessary to formulate the shooting method in such a way that its Jacobian matrix is not singular, as it might be in its most natural formulation; and it is also necessary to determine initial values for each integration which are consistent with the algebraic constraints of the DAE. These difficulties have been overcome by Lamour [16,17], who has developed a multiple-shooting code (MSHDAE) for solving fully implicit index-one differential–algebraic boundaryvalue problems, based on a multiple-shooting scheme which has a nonsingular Jacobian matrix, with dimension comparable to that used in the ODE case. His scheme also determines consistent initial values for each integration. Independently, England and Mattheij [11,22] have established the need for a dichotomically stable discretization when solving a stiff boundary-value problem in ordinary differential equations, with potentially sharp boundary layers at each end of the interval. A dichotomically stable implicit Runge– Kutta method, of order 4, has been implemented in a variable step-size initial-value integrator (SYMIRK), which could be used in a multiple-shooting approach. An explicit third-order method is also required, both as a predictor, and to provide a local error indicator, which has the correct asymptotic behaviour, both for small step-sizes in the boundary layers, and for large step-sizes where the reduced solution is required outside the boundary layers. In the case of differential–algebraic equations, of the form F(x , x, t) = 0 ∈ Rn , it is not possible to reduce the Lobatto IIIA formulae to a single equation, as can be done in the case of ODEs, but instead, a system of 2n equations must be solved for the derivatives Xi+ 1 , Xi+1 , where Xi , Xi are approximations 2
to x(ti ) and x (ti ) respectively. While the Newton iteration matrix for ODEs is of dimension n × n, in the case of DAEs it is larger, of dimension 2n × 2n, but, while the former remains nonsingular, apart from isolated values of h, for ODEs, the latter also remains nonsingular, apart from isolated values of h, for index-one DAEs. The local error indicator would then be obtained in the normal way, by comparing Xi+1 to its predicted value, but in addition to predictors for the variables Xi+1 , Xi+ 1 , it is also necessary to have predictors for 2 the derivatives. The ODE integrator (SYMIRK) has been adapted, for the solution of index-one DAEs, to produce the integrator (SYMDAE). The standard version of MSHDAE uses a BDF integrator, which is not dichotomically stable, and, as is shown in this paper, for some stiff test problems this fails to integrate across the interval of interest, while
R. England et al. / Applied Numerical Mathematics 42 (2002) 117–131
119
the dichotomically stable integrator SYMDAE encounters no difficulty. Indeed, for such problems, the modified version of MSHDAE produces an accurate solution, which does not deteriorate with increasing stiffness within limits imposed by computer word length. For some nonstiff problems, the solution is also entirely satisfactory.
2. Dichotomic stability In [5,8,10] an approach has been developed for solving stiff BVPs, which may often be considered as singularly perturbed problems with some, not necessarily explicit, small parameter ε, which typifies the thickness of possible boundary layers. For example, consider the BVP for n ODEs on the interval [a, b]: dx = x = L(t)x + q(t), dt H0 x(a) + H1 x(b) = c,
(1) (2)
where x(t), q(t), c ∈ Rn and L(t), H0 , H1 ∈ Rn×n . If some rows of the composite matrix [L(t) q(t)] are considerably larger than others, then it may be partitioned in the form: 1 12 1 1 1 11 ˜ q L L (3) [L q] = ε 21 ε 22 ε 2 q˜ L L ˜ with L(t) = O(1), q(t)
= O(1). The solution then consists of a smooth part, together with some fast modes for which x / x = O( 1ε ). In [4] it is shown that, for a well conditioned BVP (1), (2), the fundamental modes of (1) may be expressed in terms of a dichotomic basis, each member of which has either no significant growth or no significant decay. Normally the fast modes remain fast throughout the interval [a, b], and for a well conditioned problem they are negligible outside boundary layers with thickness of order ε. As with any stiff problem, it can be difficult to obtain a solution of the differential equations (1) in the large part of the interval [a, b] where the fast modes are negligible, but still control the stability of the numerical method. The situation is much more difficult than for initial-value problems (IVPs) and in [11] it is shown that it is important to use a discretization which ensures that growing modes of (1) are represented by growing discrete modes, and decaying modes by decaying discrete modes. That is the main feature of a dichotomically stable method. The following definition provides an understanding of the concept for a constant-coefficient system. Definition 1 ([11]). Let the basis solutions of the homogeneous part of an nth-order linear differential equation have components proportional to eλj t , j = 1, . . . , n. Let M h be a discretization method giving an nth-order linear difference equation with corresponding basis solutions of the discrete problem having components with growth {(rjh )i }, where rjh = 1 + hλj + O(h2 ) and i is the step number. Then M h is dichotomically stable on a region R ⊂ C if and only if, for hλj ∈ R, (1) Re(hλj ) 0 ⇒ |rjh | 1, (2) Re(hλj ) 0 ⇒ |rjh | 1.
120
R. England et al. / Applied Numerical Mathematics 42 (2002) 117–131
England and Mattheij [7] investigated methods with this numerical property, and developed [9] a code (SYMIRK) using a fourth-order symmetric implicit (Lobatto IIIA) Runge–Kutta (R–K) formula 1 1 (4) Xi+ 1 = (Xi+1 + Xi ) − h Xi+1 − Xi , 2 2 8 1 2 1 (5) Xi+1 = Xi + h Xi + Xi+ 1 + Xi+1 , 2 6 3 6 where the system of ODEs is written as x = F(x) ∈ Rn . It is important to note that these formulae can be used equally to integrate in either a forward or a backward direction, and, if satisfied exactly, will reproduce the same results in either direction. If this scheme is applied to the model problem (6) x = λ x − f (t) , |λ| 1, then Xi+1 = Xi rλh − λh
1 f 6 i
1 + 23 fi+ 1 + 16 fi+1 − 12 λh(fi+1 − fi ) 2
1 − 12 λh +
1 2 2 λh 12
(7)
where rλh =
1 2 2 1 + 12 λh + 12 λh 1 2 2 1 − 12 λh + 12 λh
is the fourth-order (2, 2)-Padé approximation to the exponential function eλh . Hence ti+1 λh Xi+1 ∼ Xi e − λf (t) dt |λh| → 0 ti
and Xi+1 ∼ Xi − fi + fi+1
|λh| → ∞ ,
(8)
so that the formulae have the correct asymptotic approximation to the solution for |λh| → 0, and, for |λh| → ∞, the numerical solution is almost parallel to the smooth particular solution, which, in the case of the model problem (6), is given to a first approximation by x(t) ∼ f (t) if |λ| 1 in a region where f (t) is smooth. Thus, any errors do not suffer from exponential growth, and the method can work efficiently using large step-sizes, in such a region. The SYMIRK code, outlined in [9], preserves the dichotomic structure of the fundamental solution space by the application of the dichotomically stable formulae (4) and (5). Not only do these formulae make it possible to integrate in either a forward or a backward direction, but they can be used to integrate accurately within the boundary layers, and also to integrate efficiently where the solution behaves smoothly if an appropriate step-size strategy is used.
3. A predictor–corrector pair To develop a variable step-size integrator based on the fourth-order, Lobatto IIIA, Runge–Kutta formulae (4), (5), another method is needed which must have order at least 3, so that the difference
R. England et al. / Applied Numerical Mathematics 42 (2002) 117–131
121
between the two values at the end of each step, which is O(h4 ), may be used as an error indicator. If such a method is to be explicit and is to avoid unstable growth over one step, then it must not use derivative values. Accordingly, an appropriate method is a four-step extrapolation formula, which, in the case of equal step-sizes is p
Xi+1 = 4Xi − 6Xi−1 + 4Xi−2 − Xi−3 .
(9)
If the 3-stage Lobatto IIIA formula is applied to the test problem (6), it produces an approximation over one step which is given by Eq. (7). Hence, for fast modes with large values of |λh|, the numerical solution is given by Eq. (8), and if this numerical approximation is substituted into the four-step extrapolator, assuming for simplicity that the four steps are of the same size, then p
Xi+1 = 4Xi − 6Xi−1 + 4Xi−2 − Xi−3
∼ Xi − fi + (4fi − 6fi−1 + 4fi−2 − fi−3 ) |λh| → ∞ (iv) + O h5 (h → 0) = Xi − fi + fi+1 − h4 fi−1
and so the principal part of the error indicator is essentially due to the interpolation error at the smooth particular solution. The SYMIRK integrator for ODEs uses this extrapolator for both prediction and error estimation.
4. Implementation and step control When solving a nonlinear system of differential equations, written in the form: x = F(x) ∈ Rn it is convenient to eliminate the intermediate value Xi+ 1 from the corrector formulae (4), (5) to produce 2 a single corrector equation:
1 1 2 h 1 (10) Xi+1 − Xi − h Xi+1 + F (Xi+1 + Xi ) − Xi+1 − Xi + Xi = 0, 6 3 2 8 6 )) is the Jacobian matrix of the differential equation, for the unknown vector Xi+1 ∈ Rn . If J(t) = ∂F(x(t ∂x then the Newton iteration matrix for this compressed corrector equation (10) is 1 1 1 (11) I − hJi+1 − hJi+ 1 + h2 Ji+ 1 Ji+1 ∈ Rn×n 2 2 6 3 12 where Ji+1 = J(ti+1 ) and Ji+ 1 = J(ti+ 1 ). 2 2 The exact value of this Newton iteration matrix will rarely be used, and so, in practice, the corrector equation will be solved by some form of modified Newton iteration, which will converge linearly, rather than quadratically as for true Newton iteration. The rate of convergence can be estimated at each iteration, by comparing consecutive norms of the residual of Eq. (10). In the SYMIRK integrator, the most pessimistic, cumulative estimate of the convergence rate is used at each iteration, and on that basis, if it appears that sufficient accuracy will be achieved within a further 3 iterations, then another iteration is taken. If it does not appear that convergence will be achieved within 3 iterations, even with an updated value of the Jacobian matrix J, then the step is abandoned and the step-size is reduced to a level at which convergence may be expected in one iteration. Convergence is considered to be achieved when the norms
122
R. England et al. / Applied Numerical Mathematics 42 (2002) 117–131
of both the residual and the Newton correction are reduced to less than 18 of the tolerance tol, which is specified as a combination of an absolute tolerance, and a tolerance relative to the current value of the solution, component by component. When convergence is achieved, the last Newton correction is not actually applied, so as to avoid the need to reevaluate the function F(x) at ti+1 and ti+ 1 . 2 Having achieved convergence to a satisfactory solution Xi+1 of Eq. (10), it is necessary to evaluate the error indicator which will be used to determine whether or not the step should be accepted, and what step-size should be used for the next step. As indicated earlier, the error indicator Ei+1 will be given as the difference between the final value Xi+1 , which satisfies Eq. (10), and the third-order predicted value p Xi+1 obtained from Eq. (9). On the first three steps of any particular solution, there will be insufficient past data to enable the predictor (9) to be used. For those three steps, it is necessary to use a lower order predictor of the same type, but no useful error indicator is then available until four steps are completed; if the fourth step is rejected, then all four steps must be rejected. Thereafter, the error indicator can be evaluated at every step, and if its norm exceeds the tolerance tol, then the step will be rejected and the step-size reduced by a factor of 0.5. If it is less than 12 tol for five consecutive steps, then the step-size will be increased by a factor
(−1/4) 1 1 if Ei+1 > tol,
Ei+1 / tol 2 32 or 2.0 if it is even smaller. If Ei+1 is between 12 tol and tol, then the step-size will be reduced by a factor
Ei+1 Ei+1
− 0.5 , 1.0 − tol tol which smoothly approximates ( Ei+1 /( 12 tol))(−1/4) for Ei+1 near 12 tol, and 0.5 for Ei+1 near tol. 5. Multiple shooting for DAEs This section will be devoted mainly to the numerical solution of DAEs of the form: F x , x, t = 0 (a t b) subject to boundary conditions G x(a), x(b) = 0,
(12)
(13)
depend only where F(y, x, t), G(xa , xb ) ∈ Rn . Let F and G be sufficiently smooth, the kernel of Fy := ∂F ∂y on t (N(t) := ker(Fy (x , x, t))), and rank Fy =: r = constant. Q(t) denotes a projection onto N(t) (see [12]). Let Q be smooth and P(t) := I − Q(t). By the assumption that the matrix A1 x , x, t := Fy x , x, t + Fx x , x, t Q(t) is nonsingular everywhere we ensure that (12) is a (uniformly) index-one DAE. Contrary to ODEs of dimension n, only r (< n) independent boundary conditions may be fixed. This is guaranteed (see [18]) iff the shooting matrix S := Gxa X(a, a) + Gxb X(b, a)
R. England et al. / Applied Numerical Mathematics 42 (2002) 117–131
satisfies the condition ker S = N(a), im S = M := im(Gxa , Gxb ),
123
(14)
where X(t, s) denotes the fundamental solution of (12). Condition (14) ensures that • (12), (13) is well-posed; • if x denotes a solution of (12), (13), then x is an isolated solution; • considering the correct description of the initial-value conditions in t = s: F x (t), x(t), t = 0, P(s) x(s) − z = 0,
(15)
and denoting the corresponding solution by x(t; s, z), then the IVP (15) has a unique solution at least in a neighbourhood P(s) x (s) − z δ of a solution x . The solution x of an IVP depends on the initial value P(s)z =: u only, i.e., x(t; s, z) = x(t; s, u). Using the subdivision of the interval [a, b]: a = t0 < t1 < · · · < tm−1 < tm = b the boundary and matching conditions, for forward multiple shooting, become G z0 , x(tm ; tm−1 , zm−1 ) = 0, Pi zi − x(ti ; ti−1 , zi−1 ) = 0, i = 1, . . . , m − 1,
(16)
with Pi := P(ti ). For every shooting point ti we need consistent initial values (yi , zi ), which are determined by the solution of Qi yi + Pi vi = 0, F(yi , ui + vi , ti ) = 0, (17) Qi ui = 0, where zi = Pi zi + Qi zi =: ui + vi (see [16,17]). A combination of (16) and (17) yields the shooting operator G(z0 , zm ) + K−1 Q0 u0 , ui+1 − Pi+1 x(ti+1 ; ti , ui ), S(ξ ) := i = 0, . . . , m − 1, Q i yi + Pi vi , F(yi , ui + vi , ti ), with ξ := (u0 , . . . , um−1 , y0 , v0 , . . . , ym−1 , vm−1 ); K must be chosen so that the image of K−1 Q0 has no intersection with Gxa .
124
R. England et al. / Applied Numerical Mathematics 42 (2002) 117–131
To solve the BVP with this shooting method we have to solve the shooting equation S(ξ ) = 0,
(18)
which has a nonsingular Jacobian with respect to ξ , in a neighbourhood of a solution. The essential part of the Jacobian matrix has the same size and also the same block cyclic structure as in the ODE case. The multiple-shooting code is designed in such a way that the directions of integration in every subinterval may be chosen independently. This algorithm is implemented in the code MSHDAE [17], which solves the IVPs (15) by means of backward differentiation formulae (BDFs) in a code modified from that of Skelboe [23]. It may readily be seen that if such an algorithm is applied to a problem such as (1), (2) in ODEs (a special case of DAEs), then the underlying BDF discretization is not dichotomically stable, and may result in an unstable numerical solution if applied to certain types of stiff problem. For DAEs in general, as discussed in [20, 21], a dichotomically stable discretization should also be used, rather than the BDF discretization used in [17], which would be excellent for stiff IVPs but not in the presence of fast growing modes.
6. A dichotomically stable integrator for DAEs The ODE integrator SYMIRK needs some modification in order to solve IVPs for DAEs. The 3-stage Lobatto IIIA method may be written as 5 1 1 X + X 1− X , (19) Xi+ 1 = Xi + h 2 24 i 3 i+ 2 24 i+1 1 2 1 (20) Xi+1 = Xi + h Xi + Xi+ 1 + Xi+1 , 2 6 3 6 and on substitution into the DAE (12) at ti+ 1 and ti+1 , this gives the following system of equations for 2 the unknown derivatives Xi+ 1 , Xi+1 : 2
5 1 1 F Xi+ 1 , Xi + hXi + hXi+ 1 − hXi+1 , ti+ 1 = 0, 2 24 3 24 2 2
1 2 1 F Xi+1 , Xi + hXi + hXi+ 1 + hXi+1 , ti+1 = 0. 6 3 6 2
(21)
It is not possible to eliminate Xi+ 1 , and so reduce this system to a single equation such as (10), and so 2
the larger system in R2n must be solved. The iteration matrix for modified Newton iteration, with respect to the derivative values Xi+ 1 , Xi+1 is 2
Fyi+ 1 + 13 hFxi+ 1 2
2 hFxi+1 3
2
1 − 24 hFxi+ 1 2
Fyi+1 +
1 hFxi+1 6
,
which is nonsingular for an index-one system of DAEs, apart from isolated values of h.
(22)
R. England et al. / Applied Numerical Mathematics 42 (2002) 117–131
125
In addition to predictors for Xi+1 , Xi+ 1 , it is now necessary to have predictors for the derivatives, 2 which, in the case of equal step-sizes, are as follows: p
Xi+1 = 4Xi − 6Xi−1 + 4Xi−2 − Xi−3 , 35 35 21 5 p Xi+ 1 = Xi − Xi−1 + Xi−2 − Xi−3 . 16 16 16 16 2 The heuristics used in the ODE integrator SYMIRK also need some slight modifications for use in the case of DAEs. The use of the Lobatto IIIA formulae (19), (20) has an advantage, over a comparable Gaussian method, that the order of consistency is the same for both the differential and the algebraic variables, and so there is no need to treat them separately. At each iteration, the formulae (5), (4) should be used to determine Xi+1 , Xi+ 1 , and the residuals for (21) should be reduced until they are smaller 2 than the local error tolerance. If the Newton step is (!Xi+ 1 , !Xi+1 ) then, in view of the form of (20), 2
the corresponding change in Xi+1 is given by h( 23 !Xi+ 1 + 16 !Xi+1 ), and so the appropriate norm of 2 this expression should also be reduced until it is smaller than the local error tolerance. The local error indicator would then be obtained in the normal way, by comparing Xi+1 to its predicted value. The ODE integrator SYMIRK has been adapted as outlined above, and the resulting integrator (SYMDAE) has been successfully used to solve a number of index-one DAE problems. The new integrator SYMDAE has also been inserted into the multiple-shooting code (MSHDAE), to facilitate the solution of differential–algebraic BVPs.
7. Test results While a number of tests have been carried out, the problem which has been studied most extensively may be described as follows:
2 t dx 1 1 2 3 − 9x + 12x + 5 x − exp 20 − x 5 − 2 = 0, dt µ
t dx 2 1 2 4 + 12x + 9x + 5 x − 2 exp −20 = 0, dt µ
dx 2 dx 3 t 1 dx 1 (23) 1 − 12 +5 − 100 3x + exp 20 = 0, 9 dt dt dt µ µ
1 2 4 dx dx t 2 dx 2 +9 −5 + 100 3x − exp −20 = 0, 12 dt dt dt µ µ x 3 x 4 − x 5 = 0, subject to boundary conditions: x 1 (0) = 0, where b = µ/20.
x 2 (b) = 0,
x 3 (b) = e,
x 4 (0) = 2,
126
R. England et al. / Applied Numerical Mathematics 42 (2002) 117–131
This problem has a known exact solution given by x 1 (t) x 2 (t) x 3 (t) x 4 (t) x 5 (t)
= = = = =
0, 0, exp(20t/µ), 2 exp(−20t/µ), 2,
(24)
but any perturbation from that solution, at one point, causes a perturbation which grows in both directions away from that point, like exp(±15t). As the stiffness parameter µ is increased, the behaviour of the exact solution is simply stretched out over the interval of increasing length, while the behaviour of the modes exp(±15t) becomes faster and faster in comparison to that of the exact solution involving exp(±20t/µ). Initial-value integrators may have problems during the multiple-shooting process, either because the initial estimates of the solution are not close enough to the exact smooth solution, or because of the accumulation of local truncation errors which are not controlled by the stability of the integrator. In either case, once there is a significant component of the fast growing modes, an integrator is going to need to take very small steps in order to represent it adequately. The idea behind the dichotomically stable integrator (SYMDAE), is that it should be able to control the accumulation of the local truncation errors and the numerical growth of the fast growing modes, so long as the component of those modes remains small, of the same order of magnitude as the local error tolerance for the integration. Results have previously been reported, in [7,9], for a related system of ODEs, using the integrator SYMIRK. Most of the tests on this system of DAEs were carried out using a single shooting interval, with integration in the direction from t = 0 to t = b = µ/20. Initial estimates of the solution were given at t = 0, and to provide a realistic test of the integrators, these were perturbed from the exact solution by an amount of the same order of magnitude as the local error tolerance, typically 10−6 . In the majority of cases, the initial estimates for [x 1 (0), x 2 (0), x 3 (0), x 4 (0)] were [10−4 , 10−4 , 1, 2], though in one case, the first two components were reduced to 10−5 . The initial estimate for x 5 (0) was given as 3, since the initial condition implies that no value of x 5 (0) is in fact used for the integration; however, an initial estimate of some kind is needed for the computation of consistent initial values. For nonstiff problems, with µ in the range from 1 to 10, the two integrators, BDF and SYMDAE, produced very similar numerical results, with final errors of the same order of magnitude (although a little larger than the local error tolerance), even when the initial estimates were perturbed by as much as 10−2 . Usually, a first integration approximated a perturbed solution, while after a single Newton correction to the vector ξ in Eq. (18), a final solution was obtained. For these cases, the dichotomically stable integrator required more function evaluations than the BDF integrator. For slightly more stiff problems, with µ in the range from 100 to 1000, both integrators ran into considerable difficulties, and were generally unable to complete an integration over the shooting interval. To resolve this problem, the interval [0, b] was subdivided into a number m of shooting intervals, as shown in the tables below, thus using a genuine multiple-shooting method. Up to a point, this produced a satisfactory solution, but the condition number of the Jacobian matrix of Eq. (18) was large in these cases, and the shooting algorithm was having increasing difficulties. For larger values of µ, the behaviour of the BDF integrator only deteriorated, as might be expected of a discretization which is not dichotomically stable in the presence of fast growing modes.
R. England et al. / Applied Numerical Mathematics 42 (2002) 117–131
127
However, for µ in the range from 104 to 1018 , excellent results were obtained using SYMDAE, the errors being well within the local error tolerance, with the integrator itself giving increasing accuracy as the stiffness µ increased. Nevertheless, for values of µ greater than approximately 1015 , the integrator SYMDAE began to encounter difficulties. The condition number of the matrix (22) is then of the order 1036 , and, owing to the finite precision of the computer, it becomes impossible to solve Eq. (21) to within the local error tolerance, without a restriction on the maximum step-size used by the integrator. This in turn causes growth of the unwanted fast modes, and makes a satisfactory solution impracticable for µ > 1018 . The principal results are given in Tables 1 and 2. In previous work on ODEs [5,8,10], it has been shown that in addition to a dichotomically stable integrator, it is necessary to develop techniques for localizing a smooth particular solution (a “pathfinder”), and for selecting appropriate shooting intervals. These aspects of the problem have not been dealt with in this paper, and no claim is made to have produced a general purpose code. Nevertheless, multiple shooting has the potential advantage over other methods, such as collocation, that it is not necessary to perform iterations to solve linear problems, such as that given here. The problem described above is a fully implicit index-one problem. In order to solve such a problem using the collocation code COLDAE [1], it is generally necessary to express it as a semi-explicit index-two problem, thereby introducing unnecessary instabilities. A case of such an instability has been Table 1 Results with the BDF integrator ν
µ = 10ν
b = µ/20
m
Iterations
Function calls
Defect
Cond(Jacobian)
0 1 2
1 10 100
0.05 0.50 5
7.7e−9 1.6e−7
2.9 2.6e+3
7720
3.9e−4
1.3e+1
1e3
50
1 3 fails 1 fails
769 1044
3
1 1 1 50 1
Table 2 Results with the dichotomically stable integrator SYMDAE ν
µ = 10ν
b = µ/20
m
Iterations
0 1 2
1 10 1e2
5e−2 0.5 5
1 1 1 10
3 4
1e3 1e4
5 10 15 17 19
1e5 1e10 1e15 1e17 1e19
1 2 fails 11 fails fails 1 1 1 3 3 fails
50 5e2 (x1,2 (0) = 1e−5) 5e3 5e8 5e13 5e15 5e17
1 1 1 1 1 1
Function calls
Defect
Cond
Total
First/last iter.
3286 6041
2032/1216 3544/1216
8.3e−10 4.0e−7
2.9e0 2.6e+3
183916
34654/1976
3.4e−6
4.8e+3
2846 2846 2846 120495 10564909 >479e6
1752/1046 1752/1046 1752/1046 85876/2136 7745524/34860
7.5e−10 2.4e−9 1.2e−8 4.6e−6 2.4e−6 1.6e−8
8.9e0 1.3e0 1.0e0 1.0e0 1.2e0 1.0e0
128
R. England et al. / Applied Numerical Mathematics 42 (2002) 117–131
observed in Ascher and Spiteri [1, Example 3], who suggest that their problem is ill-conditioned. Lamour [19] solves the same problem, using multiple shooting, without any problems of instability. Nevertheless, in order to provide some comparison, this problem has been solved using COLDAE, with a tolerance on the defect of 10−6 , and the results are shown in Table 3. Because of the special structure of this problem, it is also possible to produce an index-one formulation omitting the equation 5 for the unnecessary derivative dxdt , and results are also shown for that reduced formulation. As may be seen, the general purpose code COLDAE produces a solution using fewer function evaluations, and without the restrictions on the stiffness experienced using multiple shooting. On the other hand, as might be expected, multiple shooting, even before developing the necessary techniques for a general purpose code, clearly shows that few or no iterations are required for such a problem, while COLDAE typically requires 5 or 6 iterations. The multiple-shooting code has also been applied to a problem arising in electrical network simulation. It describes a bipolar ring oscillator which was discussed in detail in [13]. Fig. 1 shows the principle network. The current sources J1 and J2 are modelled by one transistor, J3 is constant and an induction is introduced as shown in Fig. 8 in [13]. The problem concerns the autonomous oscillation of the circuit and the results of the integration of an initial value problem up to the oscillation are used as an initial guess. Because of the unknown period of the oscillation, the dimension of the problem is increased from 15 to 16. The initial guess used at the point t = 0 is given by x(0) = [−9.81e−02, −3.24e−01, −6.03e−01, −4.51e−01, −4.49e−01, −8.81e−01, −1.64e−01, −2.41e−02, −4.98e−01, −1.50e+00, −1.01e+00, −4.10e−02, 1.26e−03, 1.26e−03, 3.06e−04, 6.50e+00] and a solution is to be determined in the interval [0, 10−8 ]. The boundary conditions consist of periodic conditions, and a phase condition which fixes the first component of x (see [18]). Shooting over a single interval, both integrators computed the solution, as shown in Table 4. The computed period is 6.50222552e−8. Fig. 2 shows the voltage at node 1 and Fig. 3 the current iss (shown as Jss in Fig. 1). The subroutines used and the complete result file are available from the authors. Table 3 Results with COLDAE ν 0 1 3 5 10 50 100
Index-2 formulation
Index-1 formulation
Final m
Iter.
Function calls
Final m
Iter.
Function calls
20 20 20 20 20 20 20
4 5 6 6 6 6 6
635 745 855 905 905 905 905
20 20 20 20 20 20 20
4 5 5 5 5 5 5
634 744 744 744 644 644 644
Table 4 Results for bipolar ring oscillator problem Integrator
tol
Iterations
Defect
Cond(Jacobian)
BDF SYMDAE
1e−6 1e−6
4 3
3.6e−6 7.5e−7
1.8e+5 1.8e+5
R. England et al. / Applied Numerical Mathematics 42 (2002) 117–131
129
Fig. 1. Bipolar ring oscillator.
Fig. 2. Voltage u1 .
This problem is also fully implicit of index one, and in addition, has nonseparated boundary conditions. In order to put it in a form which could be solved using COLDAE, it would be necessary to express it as a semi-explicit index-two problem with approximately 30 variables. Once again, this would introduce
130
R. England et al. / Applied Numerical Mathematics 42 (2002) 117–131
Fig. 3. Current iss .
instabilities. However, the version of COLDAE which the authors have is unable to solve problems with more than 20 variables, and it is not a straightforward process to increase all the relevant restrictions.
8. Conclusions The numerical solution of boundary-value problems for DAEs presents all the same problems as in the case of ODEs, together with a number of additional difficulties. In particular, using a natural formulation of the multiple-shooting process, the Jacobian matrix of the resulting set of equations is singular. This has been overcome in the code MSHDAE, by combining the computation of consistent initial values, for each integration, with the classical multiple-shooting method. Additionally, the integration method, used for solving each individual initial-value problem, needs to take account of the dichotomic structure of the fundamental modes. A dichotomically stable method for index-one DAEs has been implemented in the integrator SYMDAE, making use of earlier experience in designing such an integrator for ODEs. Numerical results are presented for an example, demonstrating that within a multiple-shooting environment, a dichotomically stable integrator should provide an efficient way of solving very stiff boundary-value DAE problems. A second example shows that it also provides entirely satisfactory results for nonstiff and nonlinear problems. The multiple-shooting code MSHDAE does not yet have all the components which would be required of a general purpose code, some of which would improve its efficiency. Nevertheless, it has been shown here that, with an appropriate integrator, it is capable of obtaining good solutions of fully implicit and nonlinear index-one problems, even with nonseparated boundary conditions. Such problems usually require more iterations using the collocation code COLDAE, and can suffer from an increase both in the number of variables and, more importantly, the index, which can be the cause of instabilities.
R. England et al. / Applied Numerical Mathematics 42 (2002) 117–131
131
References [1] U.M. Ascher, R.J. Spiteri, Collocation software for boundary value differential–algebraic equations, SIAM J. Sci. Comput. 15 (1994) 938–952. [2] K.E. Brenan, S.L. Campbell, L.R. Petzold, Numerical Solution of Initial-Value Problems in Differential–Algebraic Equations, 2nd edn., SIAM, Philadelphia, PA, 1996. [3] K.D. Clark, L.R. Petzold, Numerical solution of boundary value problems in differential–algebraic systems, SIAM J. Sci. Statist. Comput. 10 (5) (1989) 915–936. [4] F.R. de Hoog, R.M.M. Mattheij, On dichotomy and well conditioning in BVP, SIAM J. Numer. Anal. 24 (1) (1987) 89–105. [5] E. Dueñas, R. England, J. López-Estrada, Multiple shooting with dichotomically stable formulae for linear boundary-value problems, Comput. Math. Appl. 38 (9–10) (1999) 143–159. [6] R. England, R. Lamour, Dichotomic stability and multiple shooting for numerical boundary-value problems in index-one differential–algebraic equations, Z. Angew. Math. Mech. 76 (Supp. 1) (1996) 87–90. [7] R. England, R.M.M. Mattheij, Discretizations with dichotomic stability for two-point boundary value problems, in: U.M. Ascher, R.D. Russell (Eds.), Numerical Boundary Value ODEs, Progress in Scientific Computing, Vol. 5, Birkhäuser, Boston, 1985, pp. 91–106. [8] R. England, R.M.M. Mattheij, A stable sequential marching approach for two-point boundary value problems, Presented at the Workshop on Numerical Analysis and its Applications, Cuarto Coloquio de Matemáticas del Centro de Investigación y de Estudios Avanzados del Instituto Politécnico Nacional, Taxco, Mexico, 1985. [9] R. England, R.M.M. Mattheij, Sequential step control for integration of two-point boundary value problems, in: J.P. Hennart (Ed.), Numerical Analysis Proceedings, Guanajuato, Lecture Notes in Math., Vol. 1230, Springer, Berlin, 1986, pp. 221–234. [10] R. England, R.M.M. Mattheij, A sequential algorithm for solving boundary value problems with large Lipschitz constants, Rep. Appl. Numer. Anal. 87-15, Eindhoven University of Technology, Eindhoven, 1987. [11] R. England, R.M.M. Mattheij, Boundary value problems and dichotomic stability, SIAM J. Numer. Anal. 25 (5) (1988) 1037–1054. [12] E. Griepentrog, R. März, Differential–Algebraic Equations and their Numerical Treatment, Teubner, Leipzig, 1986. [13] M. Günter, U. Feldmann, CAD-based electric-circuit modeling in industry: I. Mathematical structure and index of network equations, Surveys Math. Indust. 8 (1999) 97–129. [14] E. Hairer, G. Wanner, Solving Ordinary Differential Equations II: Stiff and Differential–Algebraic Problems, Springer, New York, 1991. [15] M. Kiehl, Sensitivity analysis of ODEs and DAEs—Theory and implementation guide, Optim. Methods Software 10 (6) (1999) 803–821. [16] R. Lamour, Shooting methods for transferable DAEs—Solution of shooting equation, in: Numerical Analysis and Mathematical Modelling, Vol. 24, Banach Center Publications, Warsaw, 1990, pp. 233–240. [17] R. Lamour, A well-posed shooting method for transferable DAEs, Numer. Math. 59 (1991) 815–829. [18] R. Lamour, Oscillation in DAEs, Preprint Nr. 272, Humboldt-University of Berlin, Department of Mathematics, Berlin, 1991. [19] R. Lamour, A shooting method for fully implicit index-2 differential–algebraic equations, SIAM J. Sci. Comput. 18 (1) (1997) 94–114. [20] M. Lentini, R. März, The conditioning of boundary value problems in transferable differential–algebraic equations, SIAM J. Numer. Anal. 27 (1990) 1001–1015. [21] M. Lentini, R. März, Conditioning and dichotomy in differential algebraic equations, SIAM J. Numer. Anal. 27 (1990) 1519–1526. [22] R.M.M. Mattheij, Decoupling and stability of algorithms for boundary value problems, SIAM Rev. 27 (1) (1985) 1–44. [23] S. Skelboe, INTGR for integration of stiff systems of ordinary differential equations, Instituttet for Teleteknik, Report IT 9, Lyngby, 1977.