Accuracy of Decoupled Implicit Integration Formulas.
Dynamical systems can often be decomposed into loosely coupled subsystems. The system of ordinary differential equations (ODEs) modelling such a problem can then be partitioned corresponding to the subsystems, and the loose couplings can be exploited by special integration methods to solve the problem using a parallel computer or just solve the problem more efficiently than by standard methods.This paper presents accuracy analysis of methods for the numerical integration of stiff partitioned systems of ODEs. The discretization formulas are based on the implicit Euler formula and the second order implicit backward differentiation formula (BDF2). Each subsystem of the partitioned problem is discretized independently, and the couplings to the other subsystems are based on solution values from previous time steps. Applied this way, the discretization formulas are called decoupled.The stability properties of the decoupled implicit Euler formula are well understood. This paper presents error bounds and asymptotic error expansions to be used in controlling step size, relaxation between subsystems and the validity of the partitioning. The decoupled BDF2 formula is analyzed within the same framework.Finally, the analysis is used in the design of a decoupled numerical integration algorithm with variable step size to control the local error and adaptive selection of partitionings. Two versions of the algorithm with decoupled implicit Euler and BDF2, respectively, are used in examples where a realistic problem is solved. The examples compare the results from the decoupled implicit Euler and BDF2 formulas, and compare them with results from the corresponding classical formulas.
Introduction
. The numerical solution of stiff systems of ordinary differential
equations (ODEs) requires implicit discretization formulas. The implicit algebraic
problems are usually solved by a Newton-type iteration method which involves the
solution of systems of linear equations that often turn out to be sparse. Attempts to
parallelize a standard algorithm, as the one just outlined, often lead to disappointing
efficiency because the parallel algorithm is too fine-grained or has too large a sequential
fraction.
The waveform relaxation method [1], although not developed with parallel computation
in mind, leads to efficient parallel algorithms for the class of problems where
it works well. Multirate integration can be considered an integral feature of the wave-form
relaxation method. The relaxation part of the method is a potential source of
computational inefficiency. One relaxation iteration is fairly expensive and convergence
may be slow.
The decoupled integration methods in this paper were developed as a response
to the problems encountered in parallelizing standard integration methods and the
problems with waveform relaxation. The decoupled integration methods, like the
waveform relaxation method, exploit a partitioning of a system of ODEs into loosely
coupled subsystems. The decoupled integration methods employ only one and occasionally
two relaxation iterations. Multirate integration is not quite as natural as
for the waveform relaxation method, although it is still possible, and the parallel
implementations of decoupled integrations methods will be finer grained than parallel
waveform relaxation methods. The decoupled integration methods may be more
Department of Computer Science, University of Copenhagen, Universitetsparken 1, DK-2100
Copenhagen, Denmark, e-mail: stig@diku.dk
efficient on a sequential computer than standard integration methods, just like the
waveform relaxation method.
A previous paper introduced the decoupled implicit Euler method [2]. The existence
of a global error expansion was proved under very general choice of step size,
thus permitting the use of Richardson extrapolation. A sufficient condition for stability
of the discretization was given in [3]. This condition is called monotonic max-norm
stability, and it guarantees contractivity. Partitioned systems of ODEs are in qualitative
terms characterized as monotonically max-norm stable if each subsystem is
stable and if the couplings from one subsystem to the others are weak.
This paper is organized as follows. Section 2, "Partitioned systems of ODEs and
decoupled discretization formulas", gives preliminaries and definitions, including the
definition of monotonic max-norm stability and the presentation of the decoupled implicit
Euler and the decoupled implicit second order backward differentiation formula
(BDF2).
Section 3, "Error bounds", gives error bounds for the classical and the decoupled
implicit Euler formulas. The bounds are closely tied to the monotonic max-norm
stability condition which applies only to the Euler formulas, and the bounds are not
readily generalized to the BDF2 formula.
Section 4, "Asymptotic error formulas and error estimation", includes four subsections
treating explicit formulas, classical implicit formulas, decoupled implicit Euler,
and finally decoupled BDF2 formulas. The explicit formulas are used in error estimation
and for prediction in the decoupled formulas. The asymptotic errors of the
classical formulas are the smallest achievable for the decoupled formulas and therefore
of interest. The asymptotic errors of the decoupled formulas are given for various
modes of operation, and error estimation techniques are presented.
Section 5, "Integration algorithm", first presents general principles for decoupled
integration algorithms derived from the previous analysis. Then the details of an
implementation of the decoupled implicit Euler formula is given, followed by the
minor modifications required for replacing the Euler formula with the BDF2 formula.
The implementation employs variable step size to control the local error and adaptive
selection of partitionings among two predetermined alternatives.
These integration algorithms are used in section 6, "Examples: Chemical reaction
kinetics", where a real problem is solved using both decoupled and classical versions of
implicit Euler and BDF2. The examples show excellent performance of the decoupled
formulas.
2. Partitioned systems of ODEs and decoupled discretization formulas.
Define a system of ODEs,
continuous in Y . Stable
systems of differential equations are considered stiff when the step size of the discretization
by an explicit integration method is limited by stability of the discretization
and not by accuracy. Efficient numerical integration of stiff systems therefore
requires implicit integration methods.
Let the original problem (2.1) be partitioned as
y 0y 0y 0
y 1;0
y 2;0
necessary, the
partitioning of Y will be stated explicitly as in f r (t; y
2.1. Monotonic max-norm stability. The following stability condition introduced
in [3] plays a crucial role for the stability of decoupled implicit integration
methods.
stability.
The partitioned system (2.2) is said to be monotonically max-norm stable if there
exist norms k \Delta k r and functions a rj (t; U; V ) such that
a rj (t; U; V )ku
and the following condition holds for
the logarithmic max-norm -1 (\Delta) of the q \Theta q matrix (a rj
The condition (2.4) states that the matrix (a rj ) should be diagonally dominant
with non-positive diagonal elements. For a linear problem where A the
(a rj ) matrix can be chosen as a
Theorem 3 in [3] includes further results about monotonic max-norm stability,
including possible choices of (a rj ). Monotonic max-norm stability admits arbitrarily
stiff problems.
2.2. Decoupled implicit Euler: Stability, convergence. The decoupled
implicit Euler method is defined by the following discretization of the subsystems
by the implicit Euler formula [2]:
y q;n );
and the variables ~
y i;n are convex combinations
of values in fy i;k r. The convex combinations ~
y i;n will, in general,
depend on subsystem index r, but in order to simplify notation this dependency will
not be specified explicitly.
The method is called "decoupled" because the algebraic system resulting from
the discretization of (2.1) by Euler's implicit formula is decoupled into a number of
independent algebraic problems. The decoupled implicit Euler formula can be used
as the basis of parallel methods where (2.5) is solved independently and in parallel
for different r-values. The method can be used in multirate mode with h r;n 6= h j;n for
r 6= j, and the multirate formulation can be used in a parallel waveform relaxation
method [1].
The sequential solution of (2.5) for on a single processor will, in
general, be computationally cheaper than solving the complete system
Therefore the decoupled implicit Euler method may also be an attractive
alternative to the classical Euler formula on a sequential processor even when there
is no multirate opportunity.
Theorems 4 and 5 in [3] assure the stability of the discretization by (2.5) of a
monotonically max-norm stable problem and the convergence of waveform relaxation.
The stability condition poses no restrictions on the choice of the step sizes h r;n .
The convexity of ~ y r;n is necessary in the stability theorem of the decoupled implicit
Euler method. A convex combination would typically be a zero-order interpolation,
~
multirate organizations, where a linear interpolation is an
option.
The multirate aspect is discussed extensively in [2] and [3]. Most of the results
and techniques in this paper can readily be extended to multirate algorithms. How-
ever, presenting results for multirate algorithms would add substantially to notational
complexity, and furthermore, the algorithms and the example presented in sections
5 and 6, respectively, do not exploit multirate formulation. Therefore, the multirate
aspect will not be explored further.
2.3. Organization: Jacobi, Gauss-Seidel, relaxation. The definition of the
decoupled implicit Euler formula in (2.5) suggests a Jacobi-type organization of the
computation which is well suited for parallel computation.
Define the function G J (t; Y; ~
Y ) from the partitioning of F given in (2.2),
J;r (t; Y; ~
where
y q ). Assuming the same step
size hn for all subsystems, the decoupled implicit Euler formula can now be expressed
in the compact form,
Yn
The definition of G J given above corresponds to the Jacobi organization of the
computation. A similar GG\GammaS can be defined for the Gauss-Seidel organization,
G\GammaS;r (t; Y; ~
y q
The decoupled implicit Euler formula based on the Gauss-Seidel organization
in general, be more accurate than (2.7). The
Gauss-Seidel organization, although inherently sequential, may still be used in parallel
if the partitioning (2.2) admits a red-black reordering or a similar reordering based
on more colors (see, e.g. [4]). In the following the generic function G, meaning either
G J or GG\GammaS , will be used.
The decoupled implicit Euler formula can be used with relaxation iterations,
Y [m+1]
Yn . The computational cost of each iteration of (2.9) is approximately
the same as the cost of computing Yn from (2.7).
If the relaxation is carried to convergence, the resulting discretization is the classical
implicit Euler discretization.
When the numerical solution at t n has been computed and accepted, it is usually
used for computing the next step, and it is denoted by Yn , where Yn := Y [m]
n .
In the relaxation (2.9), each subsystem r is solved in each sweep, usually by a
Newton-type iteration carried to convergence. A relaxation iteration based on New-
ton's method with Jacobi or Gauss-Seidel organization as shown below may converge
just as fast and with substantially less computation per iteration.
@y r;n
y [m]
q;n
for
y r;n . Again the potential for parallelisation is less with
Gauss-Seidel organization than with Jacobi organization.
2.4. Decoupled BDF2. The BDF2 can also be used in a decoupled mode,
where
Stability results for this formula are only known for linear problems and constant
step size [5].
The decoupled BDF2 can, of course, be relaxed just as the Euler formula in (2.9)
or (2.10):
Y [m+1]
3. Error bounds.
3.1. Decoupled implicit Euler. The error of the decoupled implicit Euler formula
can be bounded as follows. Consider the local truncation error for subsystem r,
and the decoupled implicit Euler formula using Jacobi organization,
Y r
where ~
Y r
y q;n ).
Assume that the monotonic max-norm stability condition (2.3), (2.4) is fulfilled
with Y (t n ); ~
Y r
, and introduce the simplified notation a n
Y r
Then subtraction leads to
Y r
a n
A bound for the local error y r of the decoupled implicit
Euler formula is then (cf. Lemma 2.2, Section III.2. in [6]).
rr
a n
since
rr by (2.4) and '
Y r
3.2. Classical implicit Euler. A similar local error bound can be established
for the classical implicit Euler formula expressed as follows in a notation analogous
to (2.5):
Using the monotone max-norm stability condition, we obtain
r
where
and m1 (t) is defined as
A bound for the local error y r of the classical implicit
Euler formula analogous to (3.2) is then
are valid for all
values of step sizes but may be most interesting and useful for values where L(t n
The main difference between (3.2) and (3.3) is the last term in (3.2). The local
truncation error norm L(t n ) is O(h 2
). The order of the last term in (3.2) is O(h 2
since
If the off-diagonal elements of (a ij ) are very small, the last term of (3.1)
may be negligible, and the local error bound of the decoupled implicit Euler formula
is almost equal to the local error bound of the classical implicit Euler formula.
4. Asymptotic error formulas and error estimation.
4.1. Explicit formulas. The numerical solution of stiff systems of ODEs requires
implicit discretization formulas. A numerical integration algorithm typically
also includes an explicit formula to be used in error estimation and for computing an
initial value for the iterative (Newton-type) method used in solving the implicit discretization
problem. The decoupled implicit Euler and BDF2 formulas furthermore
require the computation of ~
Yn , where the use of an explicit formula is an option.
The explicit Euler formula Y e
is an obvious choice in
connection with the implicit Euler formula, as with
the implicit decoupled Euler formula. However, explicit formulas including F (t; Y )
should generally be avoided in connection with stiff problems when hk@F=@Y k AE 1.
The bound kY e
may be approached if Yn\Gamma1 is off
the smooth solution, which is the case if solved
approximately for Yn\Gamma1 . Therefore polynomial interpolation formulas are preferred
as predictors [7].
The linear interpolation formula,
where has the local error expansion
for
The second order polynomial predictor formula is
Y p2
where
The local error expansion for this formula is
assuming that
4.2. Classical implicit Euler and BDF2. With the classical implicit Euler
the local error Y can be expressed by
for
where all partial derivatives above and in the rest of section 4 are computed at Y (t n ).
The e-terms are obtained by substituting the expansion for Yn (4.4) into the Euler
formula (4.3), Taylor expanding Y (t at the point
finally identifying e-terms of equal power of hn .
The principal local error term h 2
can be estimated as
is the divided difference of the values Yn\Gamma2 ;
The classical BDF2 formula (cf. (2.11))
has the local error expansion
assuming that
The principal local error term of BDF2 can be estimated from
The error estimates (4.5) and (4.7), based on divided differences, are asymptotically
correct [8] if the step size varies according to
and OE(t) is sufficiently smooth. This is essentially the step size variation attempted
by the algorithms in section 5.2.
4.3. Decoupled implicit Euler. The local error of the decoupled implicit Euler
formula (2.7) depends on the definition of ~
Yn and the number of relaxation iterations
being performed (cf. (2.9)). Two different modes corresponding to different choices of
~
Yn are considered for different number of relaxation iterations.
If the partitioned system (2.2) is monotonically max-norm stable, then mode 1
discretizations are stable while the mode 2 discretizations cannot be guaranteed to be
stable.
4.3.1. Mode 1, ~
. The local error is expressed as follows:
The e [m]
r -terms of (4.8) are found analogously to the e r -terms of (4.4) from (2.7),
Y [1]
e [1]
@ ~
Y
e [1]
@ ~
Y
@ ~
The relation (@G=@Y
exploited in the expression for e [1]
3 .
Element i of the vector (@ 2 G=@ ~
evaluated as Y 0
The principal local error term h 2
can only be estimated using Y p
(cf.
)k. This inequality may be fulfilled when the
subsystems are loosely coupled, but it is by no means guaranteed by the monotonic
max-norm condition. Error estimation based on Y p
should only be used if the
quality of this estimate is somehow monitored continually.
The e [2]
r -terms of (4.8) are found from Y [2]
like the
analogous terms above:
e [2]
e [2]
@ ~
Y
The relation (@G=@Y
exploited in the expression for e [2]
3 .
For
of the classical implicit Euler formula. Therefore the
principal local error term h 2
2 can be estimated using (4.5) with Y [m]
substituted
for Yn .
The e [m]
1 and e [m]
2 terms are unchanged for further relaxation iterations, m - 3,
and
e [m]
which is identical to e 3 for the classical implicit Euler formula.
4.3.2. Mode 2, ~
. The local error is expressed as follows:
Y [m]
e [m]
e [m]
e [m]
The - e [1]
-terms for -
Y [1]
Y [1]
are identical to e 1 and e 2
(classical implicit Euler), respectively, while
e [1]
@ ~
Y
The computational cost of Y p
n is in general much less than the cost of Y [1]
n in
mode 1, and although they lead to different e 3 -expressions, - e [1]
3 , one value will
not be smaller than the other in general. However, the use of ~
n instead of
~
may compromise the stability of the decoupled implicit Euler formula; cf.
section 2.2.
As in mode 1, - e [m]
are unchanged by further relaxation iterations, m - 2
and - e [m]
2.
The development of the error terms e [m]
3 and -
e [m]
3 is given to illustrate the influence
of increasingly accurate values of Y [m]
. If the decoupling of the original problem into
subsystems is efficient, mode 1 or 2 of the decoupled implicit Euler formula gives
results very close to those of the classical implicit Euler formula and
2:
If the decoupling is poor, the differences in the higher order e [m]
r - or - e [m]
r -terms
will be significant and kY [m]
k.
The error expansions do not include multirate integration. The derivation of
error expansions analogous to the above covering multirate integration could have
been done using the basic definition (2.5), but it would be notationally complicated,
and the results are essentially the same.
4.4. Decoupled BDF2. The decoupled BDF2 (2.11) admits more "natural"
choices for ~
Yn than the decoupled implicit Euler formula. Three modes using increasingly
accurate ~
Yn are presented.
Mode 1 is expected to possess the best stability properties among the three different
modes although few theoretical results are available [5]. Mode 3 with its second
order predictor value for ~
Yn is expected to be the weakest in terms of stability properties
while mode 2 is somewhere in between.
4.4.1. Mode 1, ~
. The errors of the decoupled BDF2 are considered
for
Y [1]
and subsequent relaxation iterations (2.12). The local error is expressed as follows for
constant step size,
The e [m]
r -terms of (4.11) are found analogously to the e r -terms of (4.4):
e [1]
@ ~
Y
e [1]
@ ~
Y
@ ~
Y
@ ~
The full order of accuracy of the BDF2 is not reached so another relaxation
iteration is performed:
Y [2]
The error terms are now e [2]
e [2]
@ ~
Y
If kY (3) (t n )k AE k2(@G=@ ~
then the principal local error term can be
estimated using formula (4.7). The inequality may be fulfilled when the subsystems
are loosely coupled, but it is not guaranteed by the monotonic max-norm stability
condition.
Since the error resulting from (4.10) is O(h 2 ), this step might be replaced by the
decoupled implicit Euler formula mode 1, the result of which is denoted by Y e1[1]
n in
this subsection. The BDF2 relaxation (4.12) is then replaced by
Y [2]
Y [2]
The local error expansion for constant step size is
Y [2]
9 h 3
@ ~
Y
@ ~
Y
The decoupled Euler-BDF2 combination is expected to be the O(h 3 )-error de-coupled
formula with the best stability properties. This expectation is based on the
fact that the decoupled implicit Euler formula mode 1 is stable when the partitioned
system (2.2) is monotonically max-norm stable. Such a result does not exist for the
decoupled BDF2 formula (4.10).
Yet another relaxation iteration from Y [2]
n or -
Y [2]
Y [3]
leads to a local error expansion with the same principal local error term as the classical
BDF2 (4.6), including the case of variable step size. The principal local error term
can thus be estimated using formula (4.7).
The computational cost of using mode 1 is rather high, so therefore the following
modes are of practical interest.
4.4.2. Mode 2, ~
. Another possible choice for ~
Yn is ~
computed
from (4.1). The corresponding local error expansion for constant step size is
Y [1]
9 h 3
@ ~
Y
assuming that The principal local error term
can be estimated using formula (4.7) when the subsystems are loosely coupled.
Another relaxation iteration would lead to the same principal local error as for
the classical BDF2 so that the principal local error term can be estimated using (4.7).
4.4.3. Mode 3, ~
. Using this mode we obtain the same principal local
error term as for the classical BDF2, and therefore we also have the same possibility
of estimating the error using formula (4.7).
The use of ~
n with the decoupled BDF2 may give rise to some concern
about the stability of the discretization.
4.4.4. Summary of decoupled BDF2 formulas. The following table summarises
the local error results of the decoupled BDF2 formula. Any combination of
mode and number of relaxations having the error O(h 3
have the principal local
error C 3 h 3
after one or more additional relaxations.
~
Yn
Relax'ns
5. Integration algorithm.
5.1. General principles. The previous sections have presented the decoupled
(section 2.2) and various iteration techniques for approaching
the classical implicit Euler formula (section 2.3). Asymptotic local error expansions
have been given for different modes of employment and corresponding local error
estimation techniques (section 4.3). Finally, similar results are presented for the
decoupled BDF2 (sections 2.4 and 4.4).
All of these components can be used to construct numerous integration algorithms
where the design decisions may be guided by the properties of the problem to be
solved. The main objective of an integration algorithm is to solve a problem to a
specified accuracy using as few arithmetic operations as possible.
The following discussion only deals with the decoupled implicit Euler formula
since the stability results and error bound only apply to this formula. The BDF2
formula is expected to have analogous properties, and the implementation of the
decoupled BDF2 formula is very similar to the implementation of the decoupled Euler
formula.
5.1.1. Mode. The local error bound (3.1) shows how an accurate value of ~
Yn
reduces the influence of the partitioning on the local error. According to section 4.3,
mode 2 is to be preferred over mode 1 because of the accuracy of ~
which
can be computed at little additional cost. Although kY p
0, the reverse may be true for larger values of hn . Mode 1
should therefore be preferred if
An alternative approach for improving the accuracy of ~
Yn is relaxation (2.9).
After one relaxation iteration, ~
Yn can be considered having the value ~
Relaxation does not increase the mode, and it is attractive in this respect. However,
relaxation is computationally expensive and should only be used when it is strictly
necessary.
5.1.2. Partitioning error. The error due to the partitioning is described by
the matrix (a rj ). An aggressive partitioning with few small subsystems and otherwise
scalar equations may lead to an (a rj ) matrix with relatively large off-diagonal
elements, and it may not be diagonally dominant (2.4). According to (3.1), the error
term including kY (t n
may therefore contribute significantly.
A conservative partitioning will typically have some larger subsystems to accommodate
strong couplings and to assure numerically small off-diagonal elements in
(a rj ). The error bound (3.1) clearly shows how this may lead to a decoupled Euler
formula with essentially the same error properties as the classical formula.
The Gauss-Seidel organization (2.8) takes advantage of a non-symmetric structure
of (a rj ). Assume that the partitioned system (2.2) is reordered symmetrically in
equation number and variable number to make (a rj ) as close to lower triangular as
possible, k(a rj ) r?j k AE k(a rj ) r!j k. Then (3.1) is modified as follows,
j!r
a n
j?r
a n
and the larger lower triangular a n
rj -values are multiplied with the smaller
errors while the smaller upper triangular a n
rj -values are multiplied with the
larger errors.
5.1.3. Stability. Stability is assured by the matrix (a rj ) being diagonally dominant
(2.4) and by ~
Yn being a convex combination of previous solution values. In
mode 1, where ~
the latter condition is fulfilled. The diagonal dominance
condition may be fulfilled by a sufficiently conservative partitioning.
However, the monotonic max-norm stability condition in section 2.1 is a sufficient
condition but not a necessary condition. Therefore, mode 2 may be used without
encountering stability problems and also used when the diagonal dominance condition
(2.4) is not fulfilled. A more conservative partitioning where (2.4) is fulfilled or closer
to being fulfilled will not only improve stability but most likely also accuracy; cf. the
previous section on partitioning error. However, a conservative partitioning leads to
a more computationally expensive discretization than a more aggressive partitioning.
Relaxation iterations in mode 1 do not compromise stability, and furthermore,
the monotonic max-norm stability guarantees convergence of the process.
5.2. Implementation details. The algorithm will use the decoupled implicit
Euler formula or BDF2 with variable step size and choose between two different
partitionings: an aggressive partitioning and a conservative partitioning.
The aggressive partitioning uses the smallest subsystems possible in order to
computational cost. The conservative partitioning uses somewhat larger
subsystems in order to maintain accuracy during transient solution phases.
The decoupled implicit Euler formula is used in mode 2 ( ~
relaxation iteration (Y n := -
Y [1]
used with two relaxation iterations (Y n := Y [2]
The quality of the partitioning is monitored using (4.9). The classical Euler
solution should not be computed because of the incurred cost. In mode 1, Yn in the
partitioning criterion above is therefore replaced by Y [2]
for some oe ! 1. In mode 2, Y [1]
are replaced by -
Y [1]
Y [2]
mode 1 is used with two relaxation iterations, Y [2]
n is always available, and
the cost involved in monitoring the partitioning is negligible.
Mode 2 is used with just one relaxation iteration. Therefore the cost of a step
where the partitioning is monitored is double the cost of an ordinary mode 2 step
since -
Y [2]
n is required. When -
Y [2]
n is available in mode 2, it seems obvious to return
Y [2]
n as the result of a step, Yn := -
Y [2]
n , since -
Y [2]
n supposedly is more accurate than
Y [1]
n . However, Yn := -
Y [2]
n is only returned occasionally, with Yn := -
Y [1]
n being the
common result, and in the following example this generates oscillations in the local
error estimate while failing to improve accuracy. Therefore, Yn := -
Y [1]
n is always
returned from mode 2.
The integration algorithm can now be outlined as follows.
Initialisation, step 1:
ffl choose conservative partitioning, N
ffl no error estimation
Step n:
n from (4.1) and compute -
Y [1]
n from the decoupled backward Euler
formula using mode 2 (Y [2]
computed by mode 1 is only used when the
predicted solution Y p
n is inaccurate, as explained above)
monitor then Monitor partitioning
ffl estimate the principal local error term using (4.5), " est
ffl new step size,
" tol =" est )
(aggressive partitioning) then
The step size formula averages 1 and
" tol =" est to reduce the tendency of oscillations
in step size selection.
If the step size is decreasing, the solution may be entering a transient phase which
requires the conservative partitioning. If the aggressive partitioning is being used, the
number of steps until the next Monitor partitioning is reduced by decreasing N monitor .
The partitioning is chosen using the following conditions.
Monitor partitioning
ffl in mode 2, relax the decoupled implicit Euler formula an extra time to compute
Y [2]
n . In mode 1, use Y [1]
n in the following test
Y [1]
Y [2]
Y [1]
choose aggressive partitioning
else choose conservative partitioning
if shift from conservative to aggressive partitioning then
A switch from aggressive to conservative partitioning reflects that the aggressive partitioning
is not satisfactory. A switch in the opposite direction is tentative since
it is only known that the conservative partitioning is satisfactory so the aggressive
partitioning may also be satisfactory. Therefore the first step after a switch in this
direction is monitored.
The algorithm was developed and presented for the decoupled implicit Euler for-
mula. The modifications, beyond the obvious, to accommodate the decoupled BDF2
are small. The Euler formula is replaced by the decoupled BDF2 (2.11) in mode 3,
and the first order predictor is replaced by the second order predictor (4.2) except
the first integration step which is still taken by the Euler formula and the second step
which uses decoupled BDF2, mode 2.
The principal local error is estimated using (4.7), and a normalization by fi of the
local error estimate is introduced [6, Section III.2.], " est
=fik. The
step size selection scheme is modified slightly so that averaging is only used during
increasing step sizes:
est
The step size averaging is useful for increasing step sizes to reduce the risk of instability
while it prevents a very rapid reduction in step size at the transients.
6. Examples - Chemical reaction kinetics. The example problem is the
mathematical model of the chemical reactions included in a three-dimensional trans-
port-chemistry model of air pollution. The air pollution model is a system of partial
differential equations where each equation models transport, deposition, emission,
and chemical reactions of a pollutant. By the use of operator splitting, a number of
sub-models are obtained including the following system nonlinear ODEs:
The nonlinearities are mainly products, i.e. P i and L i are typically sums of terms
of the form, c ilm (t)Y l Ym and d il (t)Y l , respectively, for l; m 6= i.
The chemistry model is replicated for each node of the spatial discretization of the
transport part. The numerical solution of a system of 32 ODEs is not very challenging
as such, but the replication results in hundreds of thousands or millions of equations,
and a very efficient numerical solution is crucial. The problem and a selection of
solution techniques employed so far are described in [9] and [10].
The system of ODEs is very stiff with the real part of the eigenvalues of the
Jacobian along the solution ranging from 0 to \Gamma8 step sizes in the range
100 to 1000. Therefore implicit integration schemes are required, and the resulting
nonlinear algebraic problem is the main computational task involved in advancing the
numerical solution one time step. A method based on partitioning the system called
the Euler Backward Iterative method is described in [9]. It can be characterized as a
discretization by the implicit Euler formula with block Gauss-Seidel iteration for the
solution of the algebraic equations of the discretization.
An ideal partitioning would involve subsystems of size one, i.e., s
(cf. (2.2)). For this chemical reaction kinetics problem, it would be particularly
advantageous since Y i is only included in L i (t; Y ) in very few equations and never
in P i (t; Y ) (by definition). Therefore L i (t; Y )Y i is in general linear in Y i , and the
solution of a scalar equation by an implicit integration formula can be performed
without iterations.
A partitioning into all scalar equations is not viable for this problem. However,
the paper [9] identifies a total of 12 equations which should be solved in blocks of
4, 4, 2 and 2 equations. With the numbering of the chemical species used in Table
B.1 in [9], the block Gauss-Seidel iteration proceeds as follows, where the parentheses
denote the blocks of equations, f(1, 2, 3, 12), (4, 5, 19, 21), 6, 7, (8, 9), 10, 14, (15,
remaining scalar equationsg. The partitioning in [9] specified above is used as
the basis of the partitioning used for the results in [11].
In [9] the Euler Backward Iterative formula is relaxed until convergence to obtain
the equivalent of the classical implicit Euler formula. The results in [11] are obtained
from the decoupled implicit Euler formula mode 2 with one relaxation. Timing results
in [11] show good efficiency for this approach.
The aggressive partitioning and ordering used in this example for both decoupled
implicit Euler and BDF2 is described by f12, (4, 5), 20 scalar equations, (1, 2, 15, 16),
remaining scalar equationsg. The subsystems of two and four equations are enclosed
in parentheses, and the rest of the equations are being treated as scalar equations.
The conservative partitioning and ordering is specified by f12, (1, 2, 15, 16, 4, 5,
remaining scalar equationsg. Except for the block of 10 equations, the
partitioning is into scalar equations.
A partitioning is considered more aggressive than another one if it has fewer
equations appearing in blocks and/or the blocks are smaller. Concerning the two
partitionings presented here, it is obvious which partitioning is the more aggressive
since f(4, 5), (1, 2, 15, 16)g ae f(1, 2, 15, 16, 4, 5, 8, 10, 29, 30)g. The parentheses are
only retained to facilitate reference to the partitioning and ordering specifications.
The aggressive partitioning is clearly more aggressive than the partitioning presented
in [9] since f(4, 5), (1, 2, 15, 16)g ae f(1, 2, 3, 12), (4, 5, 19, 21), (8, 9), (15,
16)g, while the relation to the conservative partitioning is not obvious since one has
the larger block while the other has more equations appearing in smaller blocks.
The partitioning used in [9] is presumably based on knowledge of the chemical re-
actions, while the partitionings used here are obtained with a semiautomatic method
described in paper [12]. None of these partitionings are monotonically max-norm
stable although they come close, but despite this fact, no instability problems are en-
countered. This should not be too surprising since the monotonic max-norm stability
is a sufficient but not a necessary condition for the stability of the decoupled implicit
Euler formula.
The example used in this paper differs slightly from the example in [9], but the
differences appear in the equations being treated as scalar in all the considered partitionings

6.1. Decoupled implicit Euler. Figure 6.1 shows the errors obtained using the
classical implicit Euler formula (solid line) and the decoupled implicit Euler formula
(dash-dot line) implemented as described in section 5. The two different partitionings
mentioned above have been used adaptively. The step size is chosen by the decoupled
Global
integration
error
Time
Fig. 6.1. Integration error for the classical implicit Euler (solid line) and decoupled implicit
Euler (dash-dot line).
implicit Euler algorithm to obtain a local error estimate of 10 \Gamma3 or to a minimum
step size of 90, and the classical implicit Euler formula is applied with the same step
size selection. The discrepancy between the errors is seen to be insignificant.
A reference solution is computed using a variable step size variable order (max-
imum implementation of the backward differentiation formulas [13] with
a bound on the relative local error estimate of 10 \Gamma6 . The errors presented in the
figures are the maximum relative deviations from the reference solution measured
componentwise (the values of the components vary widely in magnitude).
The time axis is in seconds, and the initial time corresponds to 6 a.m. The
model includes the influence of the sun on some of the chemical reactions, and this
leads to very distinct transients in the solution at sunrise and sunset. The minimum
integration time step of 90 seconds is too large a step to integrate the transients
accurately, and large spikes in the global integration error are seen around 7 p.m., 5
a.m. (t=105,000) next day, and 7 p.m. (t=155,000).
The transient behaviour at sunrise and sunset can to some extent be considered
a modelling artifact. Since the large local contribution to the global error does not
influence the global error at large, it is essentially ignored by introducing the minimum
time step. The observed behaviour of the global error is not uncommon for stiff
systems of ODEs.

Figure

6.2 shows the estimated principal local error. The step size is adjusted
to keep the estimate at 10 \Gamma3 , and the algorithm is quite successful except at the
transients where the minimum step size of 90 is used.

Figure

6.3 shows the resulting step size selection (upper graph) and the selection
of partitioning (lower graph). The low value indicates the aggressive partitioning, and
the high value indicates the conservative partitioning. The values on the ordinate axis
only pertain to the step size. The total number of steps is 737, and only 182 (25%)
steps are using the conservative partitioning.
e
Estimated
local
error
Fig. 6.2. Estimated principal local error for decoupled implicit Euler.
Time
Step
size
and
partitioning
Fig. 6.3. Integration step size selected by the decoupled implicit Euler formula (upper graph)
and selection of partitioning (lower graph).
During integration using the conservative partitioning, single steps that use the
aggressive partitioning can be observed. At the first step after a switch to the aggressive
partitioning, the quality is monitored, and if it is not satisfactory, the Monitor
partitioning algorithm immediately returns to the conservative partitioning.
This example demonstrates the application of an implementation of the decoupled
Euler formula in solving a non-trivial practical problem. The computational
cost is substantially less than that for the classical implicit Euler formula, and there
is no trace of instability in the solution computed by the decoupled implementation.
e
Global
integration
error
Fig. 6.4. Integration error for the classical BDF2 (solid line) and decoupled BDF2 (dash-dot
line).
Time
Step
size
and
partitioning
Fig. 6.5. Integration step size selected by the decoupled BDF2 (upper graph) and selection of
partitioning (lower graph).
6.2. Decoupled BDF2. The numerical integration of the previous section is
repeated using the decoupled BDF2. The step size is controlled to keep the estimated
local normalized error " est at 10 \Gamma3 as before. The resulting global integration error
is shown in Figure 6.4 (dash-dot line) together with the corresponding error for the
classical BDF2 using the same step size selection (solid curve). Some deviation is
noticed, but the error of the decoupled BDF2 is comparable to the error of the classical
BDF2 formula. Comparing with Figure 6.1, it is seen that the global error of the
decoupled BDF2 formula is significantly smaller than the global error of the decoupled
implicit Euler formula. The difference originates mainly from the transients, including
the initial transient, where the minimum stepsize of 90 seconds is used.
Finally, Figure 6.5 shows step size and partitioning selection similar to Figure
6.3. The maximum step size of the decoupled BDF2 is greater than 1700 which is
three times the maximum step size of the Euler formula. The necessary number of
integration steps is 311 which is 42% of the number of steps needed by the Euler
formula. The conservative partitioning is only used for 73 steps out of 311 (23%).
There is a substantial pay-off to using the decoupled BDF2 instead of the decoupled
Euler formula, since the amount of work per step is essentially the same for the two
decoupled formulas.
The performance of the decoupled BDF2 algorithm in mode 3 is very convincing,
and there is no trace of instability, although the use of the second order polynomial
predictor for computing ~
Yn is somewhat risky from a stability point of view.

Acknowledgements

. The comments by the editor and the anonymous referees
have helped choosing asymptotically correct error estimates and improving the
presentation.



--R

The waveform relaxation method for time-domain analysis of large scale integrated circuits
Methods for parallel integration of stiff systems of ODEs
Stability of backward Euler multirate methods and convergence of waveform relaxation

A connection between the convergence properties of waveform relaxation and the A-stability of multirate integration methods
Solving Ordinary Differential Equations I
The control of order and steplength for backward differentiation methods
Estimation of errors and derivatives in ordinary differential equations

A photochemical kinetics mechanism for urban and regional computer modelling
Exploiting the natural partitioning in the numerical solution of ODE systems arising in atmospheric chemistry
Partitioning techniques and stability of decoupled integration formulas for ODEs
INTGR for the Integration of Stiff Systems of Ordinary Differential Equations
--TR
parallel numerical integration;absolute stability;multirate formulas;backward differentiation formulas;euler's implicit formula;partitioned systems