- Email: [email protected]
Stability analysis of explicit multirate methods
Mathematics and Computers 0 North-Holland Publishing
in Simulation Company
XXIII (1981)
34-50
STABILITY ANALYSIS OF EXPLICIT MULTIRATE
METHODS
W. GOMM Institut ftir Rechentechnik,
Technische
Universitiit Braunschweig,
Fed. Rep. Germany
If stiff differential systems can be divided into stiff and nonstiff subsystems, each subsystem can be integrated with different step sizes. These methods are called multirate methods, and they have been successfully used in several practical problems, especially in real-time simulations. In this paper a new way to perform the stability analysis of multirate versions of linear multistep methods is presented. This stability analysis makes use of the multirate z-transform method.
1. Introduction
can be divided into the subsystems:
In many practical problems, e.g. in real-time simulation of dynamical systems, systems of stiff differential equations have to be solved with explicit methods (because of real-time constraints, [4]). The use of explicit linear multistep methods is only possible with much work: a very short step size for integration is necessary to reach numerical stability. Solving realtime problems, the computing time of one integration step may be longer than the (real-) time step. This violates the real-time constraints. In many cases only a small part of the system of differential equations is the ‘stiff part, which requires the short step size. In these cases the system can be divided into a so-called ‘nonstiff subsystem (containing the larger number of equations) and a ‘stiff subsystem. These subsystems now can be integrated with different step sizes (‘multirates’), the stiff subsystem with a short step size to achieve stability, and the nonstiff subsystem with a longer step size. These methods with different step sizes are called multirate methods [3], split methods [9], or methods with separate step size [l].
2. The principle of multirate
wo = w(t0)
W,+I = Mw,,
t,,, T; h) .
(2.3)
Here w,+r , w,, are vectors containing all the information necessary to carry out the multistep method. The multirate method can be described with M (2.3): Ynk+i+l
=MCYnk+i,
fnk+i,
T;f),i=o,..-, k- 1 , (2.4)
( bzk+k
= M(bzk,
tnk, kT; S>
with T kT
2 integration step size for the stiff subsystem, 2 integration step size for the nonstiff subsystem, k E N = multirate factor, f:= f(y,0,t)with 0 an approximation of u. A special approximation 0 of the u-values is necessary to evaluate the stiff subsystem f at the points t&+j+l, i = 0, ._., k - 1. At these points no values of IJ are computed by the integration procedure of the nonstiff subsystem. For each n E N the equations (2.4) describe an
methods
If the following n-dimensional lem w’=h(w,t),
with dim@) + dim(u) = n, then the subsystemy’ = f@, u, t) is assumed as the stiff subsystem, Let M be a known linear multistep method applied to (2.1), with the step size T:
initial value prob-
(2.1) 34
W. Gomm /Stability ynk-k
Ynk-k+2
Ynk-k+l
I
analysis of explicit multirate methods
I T
‘.
‘Ynk
I 1
Unk-k
, .’
Ynk+2
Ynk+l
I
T
T
35
I
T “nk+k
“nk
I
I
‘Ynk+k
I kT
kT
Fig. 1. Chronological sequence of the multirate method (points of time of computation).
integration rule for a cycle of the length kT. Therefore the multirate method is regarded as a method with the step size kT (see Fig. 1). If the step size T is given to achieve stability, then the number of the functional evaluation from t,k to tnk+k 1s:
in the normal case: k * funct. eval. off+ in the case of multirate k * funct. eval. off+
k * funct. eval. ofg ,
= min(Pm
of the aims of the stability analysis
The solution of the systems of differential equations which appear by the simulation of dynamical system can be divided into transient and steady-state components. Let w’(t) = h(w, t) + r(t) )
wo = w (to)
method: 1 * funct. eval. ofg .
If the multirate method is used on a serial computer, the evaluation ofg must be spread uniformly over k time steps, in each of them one evaluation off is calculated. Under real-time constraints, this procedure requires very careful programming [4]. Otherwise this problem can be solved on a parallel computer (multiprocessor system with MIMD-structure). In this case one or more processors solve each subsystem [2]. Independent, however, of the programming of the multirate method or of the choice of the computer, it is necessary to investigate the behaviour of error and stability. The order of the local error PMULT depends on the order of the local error pM of the integration method M (2.3) and on the order pi of the approximation 0 ofv: PMULT
3. Specification
PI + 1)
be such a system r s 0 the solution nents, which will For r * 0 and t + state components, function r(t). Let
of differential equations, then for describes the transient compovanish for t + m (stable system). m the solution describes the steadythe system-answer of the forcing
w’=Xw+r,
h E C,
wo=w(to),
r=r(t)
Re(h) < 0, (3.1)
be the test equation for the stability analysis. The general form of a multistep method with 1 steps
f=OliW,_iT,&d,wk-i=O
i=O
(3.2)
applied to (3.1) produces the difference equation:
.
This result is shown in [3] and with some extensions respecting the leading error coefficients in [6]. In this article the stability problem for multirate version of linear multistep methods is described and a way is found to construct the stability polynomial for these methods.
r,_i = r(t,-i)
(3.3)
(the starting values wo, . .. . w[_~ shall be known). The z-transform of the difference equation (3.3) can be found (definition and properties of the z-trans-
W. Gomm / Stability analysis of explicit multirate methods
36
form are described in [7]):
E;f=o piZ’-’
T
‘(‘) =I;&,
so that the steady-state
(cq - XT/$)
R (2) z~-~
G(r):=&e
where W(z) and R(z) are the z-transforms
of the
sequences (wkNand(~&~N. The polynomial in the denominator of (3.4) is called the stability polynomial of the linear multistep method (3.2)
g
(ai - ATPi)Z‘-j =:
p(z) - XTu(z)
.
(3.5)
One of the roots ai, j = 1, .... 1 of the polynomial (3.5) is called the principal roar (here ar), the others are the extraneous roots. It is known that ul is an approximation of eATand that T-’ In ai is an approximation OfX. The method (3.2) is called absolutely stable for AT, if lail < 1 for j = 1, .... 1. \aiJ= 1 is only possible for simple roots aj. The region of absolute stability is the set of all XT, for which the method (3.2) is absolutely stable. The method (3.2) is called rehztively stable for XT, if]a,]>ifzil,j#l. The analysis of the steady-state components of the numerical solution of the differential equation is done only with one of the most important forcing functions: r(f) = eiwf. Using the Laplace-transform the transfer function of the system (3.2) is computed
=-+(s) =-&-&.
W(s)
(3.6)
The factor
A=:
=
is
.
The factor C(iw) = (io - A)-l is called the frequency response, it represents the change of phase and amplitude of the forcing function eiwr, caused by the system. The transfer function of the numerical system is derived from (3.4) (see [5]), it is (with (3.5)):
C(z) =
T@) p(z) - hTu(z) ,=,iwT
so that the steady-state solution is
(3.8)
’
component
of the numerical
Now four aims of a stability analysis can be formulated: (1) to find the region of absolute stability, (2) to investigate the quality of the approximation of ehT by al resp. Xby T-l lnai, (3) to investigate the relative stability, (4) to investigate the quality of approximation of the exact transfer function G(L) by the numerical transfer function C(eiW3. Now it remains to find four analogous points for the stability analysis of multirate methods.
4. Stability analysis for multirate versions of linear multistep methods It is not possible to use the (one dimensional) test evation (3.1) to perform the stability analysis of multirate methods. Therefore at least a two-dimensional test equation is needed:
G(s)
is called the Luplace transfer function of the system. The output function w(t) of (3.6) is w(t)
iwt
component
&
(eiwr - eht)
+&
with the forcing functions with
8
= r(‘)(t), i = 1,2 and
eiWf
(p) let nl, nZ be the eigenvalues of for t+m,
(3.7)
Re(qi) < 0, i = 1,2
LE
with h,-l
(4.2)
W. Comm /Stability
(y) let the first equation
be the stiff equation:
Re(qi) << Re(q2) < 0. The application of the multirate version (2.4) of the general form of a linear multistep method (3.2) to the test system (4.1) leads to the difference equations:
5 j=0
ffjynk+i-i
-
h,T,$
@Ynk+i-j
-
/JT!$
T2
be two different transformation variables for k E N. The z-transform is applied to sequences with the step size kT, the zk-transform to sequences with the step size T. The following theorem shows the transfer from zk- to z-transform: Theorem 4.1 ([7,8]). If the zk-transform of Cjn)na is known:
fljhk+i-j
1 =
Y(zk) =
Pj$A+i-i
c yn zkn ,
with Cj= e2tilk
where the values fin&&j denote the approximations of u,k+j_-j, which are not computed in the second difference equation of (4.3). The difference equations (4.3) define various sequences with different step sizes: the sequence (j,JnE~ with the step size T and its subsequence with the step size kT, the sequence (unk)nCZJJ tink)nEN with the step size kT and a sequence (&k)n~ of approximated values, which has the subsequence (v,&jna. The results of the multirate method (as a method of the step size kT) are two sequences (unk)nEN of the step size kT to approximate tink)ntN> the exact values v(t), v(t) at the points (t,&N. In Section 3 the z-transform method is used to construct the stability polynomial of a linear multistep method applied to the scalar test equation (3.1). To construct the stability polynomial of the multirate method, so-called multirate z-transforms are needed, which are z-transforms with different transformation variables. The problem of these multirate z-transforms is the combination between transforms with different variables. Let
eTs,
then the z-transform of the subsequence
y(z) =; ;ij
(jnk)nEN is
Y(z&‘) (kth root of unity).
The second type of transfer from z- to zk-transform is not unique: it depends on the chosen approximation. The principle of this transfer can be described: let G(Z) = CrZo g&z-” be a known z-transform of the sequence (g&&m. Using suitable approximation methods like inter- or extrapolation polynomials, a piecewise continuous functiong*(t), t Z 0 is constructed, so that
g*(t) continuous on [nkT, (n + 1) kT) for all n E N. Using the Laplace transform method, g*(t) is transformed to L[g,(t)]
= C,(s)
= H(s) G(ekTS) = H(s) G(z)
(see Section 5 for more details). The Laplace transform H(s) is often called interpolator or extrapolator or simply called hold. The Laplace transform G,(S) is transformed to a zk-transform using known tables [7] and the following Theorem 4.2 ([6]). Let the Laplace transform Y(s) = G(z) F(s),
z = ekTs
(-z-transform) be known. If the zk-transform F(zc) of the Luplace transform F(s) can be computed, then
and zk = e Ts
zk =
n=o
for i = 1, .. .. k (4.3)
z = ekTs
31
analysis of explicit multirate methods
= Zl/k
(*
zk-transform)
y(zk)
= G(z)
F(Zk)
>
fk = eTs = zllk
W. Gomm / Srabiliry analysis of explicit multirate methods
38
If the multirate z-transform method is applied to the difference equations (4.3), then an approximation 0 of the u-values has to be found. These approximated values can be computed with extrapolation polynomials from the u- and u’-values at the points t,k-jk, j = 0, 1, ... . In the case of the linear test system (4.1) the following zk-transform can be found, formally expressed: p(zk) = Hh,(zk)
v(z) + H&k)
ycz)
.
(4.4)
In Section 5 the structure and the way of construction of these holds HA,(zk), H,(zk) Will be shown exactly. The following theorem shows the principle of construction of the stability polynomial of the multirate versions of linear multistep methbds.
form to z-transform is only of theoretical importance because the computing effort (c number of complex functional evaluations Y(z&j), j = 0, 1, .... k - 1) becomes larger for larger k. A computing rule with a computing effort independant of k has to be found. 5.1. Construction
of holds
Because of real-time constraints the holds have to fulfil the following conditions: (a) For the construction of holds, only u- and u’values can be used which belong to the past and which are used for the integration of the nonstiff equation. (b) Let t,& be the last used point of time, then the value ukk = g@nk, U&, tnk) cannot be used for the approximation of the values 1, .... k - 1, because the functional evaluation of g needs the whole interval (tnk, tnk+k). Condition (a) allows only the use of extrapolation polynomials for the construction of holds. This can be done by applying the Taylor series, so that the order of error is I- 1: fi,k+j,
Theorem 4.3 ([6]). The stability polynomial of the multirate version (4.3) of a linear multistep method (3.2), applied to the test system (4.1) is: 0 = P(z)
- PTA,(Z)1 [P(Z) - h,kTo(z)l
- /LekTZ u(z) A*,(Z)
(4.5)
ir(nkT+ T) - u(nkT+
with
j
=
T) = O(T~)
(5.1)
for 7 E [0, k7’) and n EN. The set-up for the hold (4.6)
ml O(nkT+ 7) = c
m=O
and
‘5 k B(z):=A,(z)
1
j=O
.
@ktj) P(zkt’)
He(zktj) -
A, T@kf’>
a,(T)
U,k-mk
m2 ,
(4.7)
where 5 = e2nilk is the kth root of unity and p(z), u(z) are from (3.6). k E N is the multirate factor, and the zk-transform of the approximated function 0 is given by (4.4).
5. Implementation of the stability analysis of the multirate versions of Adams-Bashforth methods The application of the Theorem 4.3 to construct a stability polynomial for a concrete multirate method needs some additional considerations: (1) Construction of holds (4.4). (2) Theorem 4.1 about the transfer of zk-trans-
+ T
c m=l
Pm(T)
ULk-mk
(54
leads to the linear system ml
( l&f
.
“1 = gl
(5.3)
[a,(7)(-mIil
+ .EI [Pm(T)(-mP1
forj=l,...,l-1. The approximation i)(t) for t > 0 can be described with the following transformation of variables t=nkT+rer=t-nkT
fornEN:
-J
W.
c c a,?(t n=o m=o
- nkT) u,k-_mk H,(s) = 0 .
m2 +
c m=1
flm(t
.
-
nkT)
t&k-mk
1
[u(t - nkT) - u(t - (n + 1)
(5.4)
Example 5.2. Hold of the order 1 = 1, group B, the set-up (see (5.2)) O(nkT+ 7) = (110(T)u,k + T& (7) t&-k
WI
results in the linear system
with t-co, t>O.
0,
u(t) = { 1,
1 = oo(r) 1
(Let the values u-i, v!_i, j = 1, .. .. max(mr, m2) - 1 be zero because the problem of starting values is not important for the stability behaviour.) The terms am(t - nkT), P,(t - nkT) are polynomials in (t - nkT). Eq. (5.4) is transformed, using Laplace transforms and the linear test equation (4.1)
*
&=;al(T)
P(s) = I&,(s)
V(ekTS) + H,(s) Y(ekTS) .
o(t) = c
n=o
[u,k + (t - nkT) &k-k] . [u(t-nkg--(t
(5.5)
Extrapolation polynomials to the degree 3 are computed in [6], they are divided into two groups: Group A: m2 = 0 in eq. (5.2) Group B: ml = 0 in eq. (5.2). It is H,(s) = 0 for all holds of the group A. Now two examples of holds are computed.
i3(nkT + 7) = 0(,,(r) u,k + @l(T) u,k-k
P(s)= I]
--t--de i
(7)
>
J&Js)
5
=f?+
H,(s) = -St
&=-al(T)
+y)vnk-Funk-k]
. [u(t - nkT) - u(t - (n + 1) k7’)] Application in
of the Laplace transform
kTs 2
)
s2 ekTs
*
Y(ekTs)
p(
_kT+z u( Sz]j sz2
with the solution (ILL = 1 + r/kT, cxl (T) = -r/kT. This solution is set into (5.4):
n=O ((1
ekrs _ 1
kT
so that with z = ekTs:
results in the linear system
1 = (110(T) •t a1
-(n+l)kT)].
Application of the Laplace transform method under consideration of u’(t) = ey(t) + &u(t) results in
.V(ekTS)tE Example 5.1. Hold of the order 1= 1, group A, the set-up (see (5.2))
&(q=j.
The solution is set into (5.4):
161:
i)(t) =
39
so that with z = ekTs:
Ml
m
i)(t) =
Gomm / Stability analysis of explicit multirate methods
method results
z-l __ s2z
’
5.2. Practical execution
of the transfer between zk-
and z-transform
Let g(t) be the time function belonging to the’Laplace transform G(s), then Zk [G(s)] is called the zktransform belonging to the sequence (g(nT)),a. Now the transfer between zk- and z-transform which is used in Theorem 4.3 can be described: The z-transform G(z) of the zk-transform G(zk) has to be found (written: G(zk) o--c G(z)) @k) P@k)
Zk [HAu(s)l -
A, T’J(Zk)
.
_ ‘h”(4 B(z)
W. Gomm /Stability
40
analysis of explicit multirate methods
-=kT
and
zkT
1)2 -m.
(zk-
It follows that The holds in Examples 5 .l and 5.2 and in [6] have a special structure, so that only some fundamental subproblems have to be solved: The calculation of the z-transforms of @k) P(zk)
Zk -
b-j]
for j = 1, 2,....
Tz(~21(4+~2o(~n)) &zkb-21
M
@k)
02r(a,)
(5.6)
’
-a,)-’
+ (1 -a&)
k-‘(1
+ (ah - 1) k-‘(1
-a;)-’
= (1 -a&)-’
,...,
!??=I ,..., 1. (5.7)
The third step requires a lot of considerations will be demonstrated by an example.
which
Example 5.3. Computation
of
&zk[s-21
of the z-transform
m = 1, .... I
7
(see also
~Sj(a,,z)
+----&zk[s-il
the z-transform
Zk [s-‘1 = Tzk(zk - 1)-’ *
1)2 =ZkT
1
1 + (1 -
(1 -a,)2(zk-
&)(Zk
-
1)2+
(1 -
a,)2(zk
- a,)
The following transfers are true [6,7] : zk -•4zk-1
zk -----e----tzk-am
[6])
1
~
Zk ~&&)I
and
Zk-am
-zk[&(~)l Zk-am
can be computed. (2) The decomposition (5.6), (5.7) into partial fractions is remultiplied with the z-transforms mentioned in (1). So the terms (4.6) (4.7) can be computed. (3) Now the polynomial (4.5) can be computed directly. of Example 5.1). The
(z - II2
TZk
@k) =
-1
Tz/c *(zk-a&(zk-
from
of
Example 5.4 (continuation zk-transform of P(s) is
gives [7]:
.
The stability polynomial (4.5) is constructed the separate ‘units’: (1) Using the zk-transforms
1 j=1,2
-u$)-~
and
where 1 is the degree of the denominator (c numbers of steps of the Adams-Bashforth methods). (3) Computation of the z-transform of the terms &zk[s-‘1,
a&)
=: S2(am, z)
(~~~(a~) = (-&)(l
This can be done in 3 steps: (1) Computation of the zk-transform Zk [s-j]. (2) Decomposition of the rational zk-transform into partial fractions
= ii% zk - a,
-
with
.
A, T(-+k)
p(zk) - A,, Tu(zk)
(z - 1)2(z
zk-
1
kT(zk-
1)2 1
v(z)
Z2
1)
1 ’
=: &,(zk)
v(z) .
WithoIl = (1 - ak,)(l -a,)-’ and QZO(U~), a2 1(am) from Example 5.2 the z-transform &
H*u(zk) O---X
[Z(~ll(%)
+ ~21(6n))
Z z-l’
-
Z z-am
k'
(a20(4
-
~ll(4?l))l
Hz
-
dwl
*
Then the stability polynomial of the multirate version of the Adams-Bashforth-2-step method with
41
W. Gomm / Stability analysis of eqolicit multirate methods
the first-order hold of group A is: G22(z)
0 = B(z) [p(z) - x”kTu(z)]
- &W%(z)
=
A (z) with (see also Theorem 4.3)
with
Det(z) = P(z)
A(z)
a!&)
(z -
B(z) = z . *i,
= ?$l bm
- N&(z)1
b(z)
-
Wb~(z)l
- wkT2u(z)A,(Z) >
[z@ll@m)+
Cl(z) =‘(z)
a21(4)
Tmel eiwim -a,
’
C,(z) = kTu(z)
with a,,,, b, from the partial fractions of (5.6). -(a20(~m)-all(am))l
fi
(z-u,3
j#”
where
6 (iw) := [r7, - G(tn)]/rZ(tn)
u(zk) m=l
zk
-am
dzk)
-
In the case of the scalar test equation (see (3.8)), the relative error between the exact transfer function and its numerical approximation is computed by
A, T"czk)
There exists an analogous equation system (5.8)
and
a(&) @k)=
,@k)=&-k,
-l&k+
:= G(eiWkT)[G(io)]-’
- 1.
= G(eiWT) [G(io)]-’ for the transfer
-E,
0.5 .
(5.9)
where
Beside the computation of the stability polynomial the transfer function of the numerical system analogous to (2.8) has to be found. It is a (2 X 2)-system of transfer functions. Using the forcing functions
[(-(ia)]-’
=
i0 - A, -’ --f
iw-
X,
1’
5.3. Numerical results
r(‘)(t) = rc2)(t) = eiWt the transfer matrix can be described by Theorem 5.1 ([6]).
The transfer matrix
1
(5.8)
has the following elements
Gl=& q*(z)
=-
1
Det (z)
Gzl (z) = __
1
Det (z)
[P(Z)-
k~“@(Z)l
Cl(Z)
The numerical stability analysis of multirate versions of Adams-Bashforth methods shows a direct connection between the size of the coupling between y and u (=^the elements p, E) and the quality of stability properties of the multirate methods. Therefore a quantitative measurement of the coupling is necessary to describe that phenomena quantitatively. Eq. (4.5) and the computation of the holds show that only the product tie is used in the equations. The size of 1-16 has to be measured in relation to the stiffness ratio, which is given by
9 8B := IRe(rlr)I/IRe(~2)l
PW,(Z)
ekTo(z)
c2 (z) 9
Cl(Z)
7
in the case of the linear equations (4.1) (4.2). The matrices A and a - A, a E R, a > 0, however, must have the same degree of coupling because these matrices have the same stiffness ratio and with the step size T or T/a the same stability polynomial. Therefore the following definition is made:
42
W. Gomm / Stability analysis of explicit multirate methods
Definition 5.2. The degree of coupling GK is measured in relation to the stiffness ratio SR GK := ,t: ,i”,“, %
(sign GK = sign&e)) .
(5.10)
The influence of the size of coupling on the stability properties can be shown in the following way: the multirate method is applied to a set of similar matrices, which are used in (4.3). These similar matrices have the same eigenvalues ql, Q, but different degrees of coupling. The comparative value is the matrix with /_E= 0, that is a degree of coupling GK = 0. In this case the multirate method changes into the application of the linear multistep method to the equation y’ = Xvv = 7,~ with the step size T and to the equation V’= X,LJ= QU with the step size kT. It is well known that the exact solution of the system (4.1) is a linear combination of eqlt and eqzr. But the component e ‘?lr is not a significant part of the solution if the first equation of (4.1) is the stiff equation (see (4.2) 7) and if t > 0. So the second investigation aim in Section 2 is only valid for the principal eigenvalue Q, which is approximated by (kT)-’ In al, where a1 is the principal root of the multirate method. Thus the relative error of the principal eigenvalue can be measured by
The numbers at the curves denote the degrees of coupling which are computed by (5.10). These numerical results show that the multirate versions of the Adams-Bashforth-2-step method have for weakly coupled systems (GK < 5%) nearly the same stability properties like the multirate method applied to the comparative matrix (GK = 0). This is also true for all degrees of coupling, but with step sizes kT < 0.1. To integrate the whole system with the normal version, however, a step size less than 0.02 is necessary to achieve numerical stability. Instead of the error matrix (5.9) the transfer functions of the both componentsy, u and their relative errors are computed to present the numerical results in a clearer way: using the forcing functions r(‘)(t) = rW(t) = eiWr the exact transfer functions can be computed (iw - h,) + P y”(t) =
U(t) =
eiwt =: G,(iw)
eimr ,
eiwf =: G,(iw)
eiwr
(iw - h,)(iw - hY) - PE (iw - X,) + E (io - X,)(io - h,) - ye
and the numerical
transfer functions
(5.11) Some numerical results for the multirate versions of the Adams-Bashforth-2-step method with first order holds and the multirate factor k = 50 are shown in the following figures. One of the hold is out of the group A, the other out of the group B (see the text in the figures). The absolute values of the principal root al (the approximation of e km2) and the absolute value of the absolute greatest extraneous root are shown in the Figs. 2-5(a). The relative errors of the approximation of q2 (computed by (5.11)) are shown in the Figs. 2-5(b). The regions of absolute stability are shown in the Figs. 6,7 (only the upper part of the symmetric regions). In all figures the symbols ‘+’ denote the values, which results from the application of the multirate method to the comparative matrix with GK = PC = 0.
The relative errors can be computed
by
IA&o)1
:= I~Y(eiWk~[Gy(i~)]-l
- 1 I,
IA,(iw)I
:= IG,(e iwkT)[G,(iw)J-’
- 1I .
Table 1 shows the errors of the following example: the multirate versions of the Adams-Bashforth-2-step method with first-order holds, the multirate factor k = 50, and the step size kT = 0.1. The results of this example show that the maximal error of the two components of the normal version differs only a little from the maximal error of the multirate version with the hold of the group A. Both, however, are better than the multirate version with the hold of the group B.
W. Gomm / Stability analysis of explicit multirare methods
43
1.0
0.8
0.6
0.4
0.2
a 0.0
1 0.0
0.1
0.2
0.3
0.4
0.5
0.7
0.6
AB- 2
OROER OF EXTRAP. A 1~ PRINCIPAL ROOT A. MAX. EXTRANEOUROOT Z-PLANE, EIGENtALUESs ( -0.500000E02 , 0.000000E00 ) 1 -O.lOOOOOEOl , o.boooooEoo
0.30 t 0.25
t
AB- 2 ORDER OF EXTRAP. A 1, REL. ERROR OF PRINCIPAL ROOT EIGENhLUES, f -0.500000E02 , 0.OOOOOOEOO ) 1 -0.100000E01
Fig. 2.
,
0.000000E00
I
I
I
0.8
W. Gomm / Stability analysis of explicit multirate methods
44
1.0
0.8
0.6
0.4
0.2
0.0
J 0.0
I
0. I
0.2
0.3
0.4
0.5
0.6
0.7
0.8
AB- 2 ORDER OF EXTRAP. A I, PRINCIPAL ROOT A. MAX. EXTRANEOUSROOT Z-PLANE, EIGENtALUES, ( -0.500000E02 , 0.000000E00 ) ( -O.lOOOOOEOl , O.bOOOOOEOO )
0.30
0.25
0.20
0.15
0.10
0.05
0.00
0.0
0.1
0.2
0.3
0.4
AE- 2 OROER OF EXTRAP. A I, REL. ERROR OF PRINCIPAL ROOT ( -O.lOOOOOEOl EIGENbALUES, I -0.500000E02 , 0.000000E00 )
Fig. 3.
0.6
0.5 ,
0.000000E00
0.7 1
0.8
W. Gomm /Stability
45
analysis of explicit multirate methods
1.0
0.8
0.6
0.4
0.2
0.0
0
)
0.1
0.2
0.3
0.4
0.5
0.6
0.7
AE- 2 ORDER OF EXTRAP. B 1, PRINCIPAL ROOT A. MAX. EXTRANEOUIROOT Z-PLANE, ( -0.100000E01 , 0. boooooEoo EIGENtALUES, ( -0.500000E02 , 0.000000E00 )
)
0.25
AS- 2 ORDER OF EXTRAP. B I, REL. ERROR OF PRINCIPAL ROOT EItENtALUESn I -0.500000E02 , 0.000000E00 I ( -O.lOOOOOEOl
Fig. 4.
,
0.000000f00
)
0.8
W. Gomm /Stability
46
analysis of explicit multirate methods
1.0
0.8 0.4
I 0.6
20
0.4
0.2
a 0.00.0
0.1
0.2
0.3
0.4
0.5
0.6
A&
2 ORDER OF EXTRAP. 6 I, PRINCIPAL ROOT A. MAX. EXTRANEOUIROOT Z-PLANE, EIGEN(rALUES, I -0.500000E02 , 0.000000E00 ) I -0.l00000E0l , o.booooofoo
0.7
0.8
)
0.30
0.25
0.20
0.15
0.10
0.05
0.00
0.0
0.1
0.2
0.3
0.4
0.5
AS- 2 OROER OF EXTRAP. S 1, REL. ERROR OF PRINCIPAL ROOT ( -O.lOOOOOEOl EIGENtALUES, I -0.500000E02 , 0.000000E00 )
F3g. 5.
0.6 ,
0.000000E00
0.7 )
0.8
W. Gomm / Stability analysis of explicit multirate methods
41 1.0
a
AB- 2,
ORDER
OF
EXTRAP.
A
I,
REGION
OF STAB.,
STIFFNESS
RATIO
50
1.c
AEI- 2, ORDER
OF
EXTRAP.
A
1,
REGION
OF STAB.,
Fig. 6.
STIFFNESS
RATIO
50
W. Gomm / Stability analysis of explicit multirate methods
48
1.0
-1 AB-
:o
2, ORDER
OF EXTRAP.
Bl,
REGION
OF STAB.,
STIFFNESS
RATIO
50
1.0
b
AB-
2, ORDER
OF EXTRAP.
01,
REGION
OF STAB.,
Fig. I.
STIFFNESS
RATIO
50
W. Gomm / Stability analysis of explicit multirate methods
Table 1 Adams-Bashforth-2-step GK
+I%
ilO%
?50%
w
method,
k = 50, kT=
Normal
49
0.1
version
1. order hold B
1. order hold A
IA,(iw) I
I A,(iw)
0.0 0.0
0.33E-1 0.33E-1
0.80E-3 0.82E-3
0.85E-3 0.83E-3
0.47E-3 0.47E-3
0.35E-1 0.35E-1
0.73E-3 0.1 lE-2
0.29E-2 0.29E-2
0.35E-2 0.34E-2
0.29E-2 0.29E-2
0.37E-1 0.37E-1
0.25E-2 0.35E-2
0.25E-5 0.24E-5
0.35E-1 0.35E-1
0.21E-1 0.21E-1
0.35E-1 0.35E-1
0.43E-1 0.43E-1
0.35E-1 0.36E-1
0.33E-5 0.32E-5
O.lOEO O.lOEO
0.41E-1 0.40E-1
O.lOEO O.lOEO
0.50E-1 0.50E-1
O.lOEO O.lOEO 0.77E-2 0.84E-2
IA,(b) I
Ih&w) I
IA&w)
0.0
0.0 0.0
0.0 0.0
0.0 0.0
0.5
0.77E-7 0.75E-7
0.47E-3 0.47E-3
1.0
0.43E-6 0.43E-6
3.0
5.0
I
I A,(k)
I
-
0.0
0.0 0.0
0.0 0.0
0.0 0.0
0.0 0.0
0.60E-1 0.52E-1
0.5
O.l3E-6 O.l3E-6
0.47E-3 0.47E-3
O.l5E-2 O.l3E-2
0.49E-3 0.50E-3
0.62E-1 0.5 7E-1
0.73E-2 0.84E-2
1.0
0.81E-6 0.78E-6
0.29E-2 0.29E-2
0.67E-2 0.62E-2
0.30E-2 0.30E-2
0.69E-1 0.67E-1
0.5 7E-2 0.1 OE-1
3.0
0.72E-5 0.67E-5
0.35E-1 0.35E-1
0.5 3E-1 0.50E-1
0.35E-1 0.36E-1
O.llEO 0.89E-1
0.27E-1 0.43E-1
5.0
0.14E-4 O.l2E-4
O.lOEO O.lOEO
0.12EO O.llEO
O.lOEO O.lOEO
0.14EO 0.14EO
0.94E-1 O.llEO
0.0
0.0 0.0
0.0 0.0
0.0 0.0
0.0 0.0
0.94E-1 0.36E-1
0.37E-1 0.44E-1
0.5
O.l6E-6 O.l5E-6
0.47E-3 0.466-3
0.24E-2 O.l3E-2
0.94E-3 0.97E-3
0.95E-1 0.50E-1
0.36E-1 0.43E-1
1.0
O.lOE-5 0.97E-6
0.29E-2 0.29E-2
0.94E-2 0.70E-2
0.38E-2 0.41E-2
O.lOEO 0.73E-1
0.32E-1 0.43E-1
3.0
0.1 lE-4 O.lOE-4
0.36E-1 0.35E-1
0.80E-1 0.74E-1
0.35E-1 0.39E-1
0.15EO OS7EO
O.llE-1 0.76E-1
5.0
0.26E3-4 0.23E-4
O.lOEO O.lOEO
0.20EO 0.18EO
0.97E-1 O.llEO
0.23EO 0.25EO
0.63E-1 0.15EO
6. Discussion The necessity for a two-dimensional test equation is the essential point of difference between this analysis and the well-known stability analysis of the ‘normal’ multistep methods. We have shown a general construction rule for the stability polynomial of the multirate versions of linear multistep methods, using the
I
multirate z-transform method. This rule has been applied to the stability analysis of multirate version of the Adams-Bashforth methods and has been demonstrated in detail on the Adams-Bashforth-2-step method. The stability properties of these multirate versions differ for weakly coupled systems slightly from the stability properties of the methods applied to the comparative matrices.
50
W. Gomm /Stability
analysis of explicit multirate methods
References [l]
H. Blum, Numerical integration of large scale system, using separate step size, in: A.W. Bennett and R. Vichnevetsky, Eds., Numerical Methods for Differential Equations and Simulation (North-Holland, Amsterdam, 1978) 19-23. [2] W. Bub, Private communication. [ 3) C.W. Gear, Multirate methods for ordinary differential equations, Department of Computer Sciences, University of Illinois, Report UIUCDCS-F-74-880. [4] C.W. Gear, Simulation: conflicts between realtime and software in: J.R. Rice, Ed., Mathematical Software III (Academic Press, New York, 1977) 121-138.
[ 5I E.G. Gilbert, Dynamic-error
analysis of digital and combined analog-digital computer systems. Simulation 6 (1966)241-257. [6] W. Gomm, Stabilit&suntersuchung van expliziten Multirate-Methoden zur numerischen Liisung von Anfangswertproblemen bei gewljhnlichen Differentialgleichungen, Dissertation, TU Braunschweig (1980). of the z-Transform (7 E.I. Jury, Theory and Application Method (Wiley, New York, 1964). ]8 I M.A. Neighbors and CA. Halijak, Effective computation of exact multirate z-transforms, in: Vogt and Mickle, Eds., Modelling and Simulation Vol. 8, Past 1, Proc. 8th Annual Pittsburg Conference (1977) 539-542. methods for simultaneous [9 J.R. Rice, Split Runge-Kutta equations, J. Res. Nat. Bureau of Stand. 64B (1960) 151-170.
in Simulation Company
XXIII (1981)
34-50
STABILITY ANALYSIS OF EXPLICIT MULTIRATE
METHODS
W. GOMM Institut ftir Rechentechnik,
Technische
Universitiit Braunschweig,
Fed. Rep. Germany
If stiff differential systems can be divided into stiff and nonstiff subsystems, each subsystem can be integrated with different step sizes. These methods are called multirate methods, and they have been successfully used in several practical problems, especially in real-time simulations. In this paper a new way to perform the stability analysis of multirate versions of linear multistep methods is presented. This stability analysis makes use of the multirate z-transform method.
1. Introduction
can be divided into the subsystems:
In many practical problems, e.g. in real-time simulation of dynamical systems, systems of stiff differential equations have to be solved with explicit methods (because of real-time constraints, [4]). The use of explicit linear multistep methods is only possible with much work: a very short step size for integration is necessary to reach numerical stability. Solving realtime problems, the computing time of one integration step may be longer than the (real-) time step. This violates the real-time constraints. In many cases only a small part of the system of differential equations is the ‘stiff part, which requires the short step size. In these cases the system can be divided into a so-called ‘nonstiff subsystem (containing the larger number of equations) and a ‘stiff subsystem. These subsystems now can be integrated with different step sizes (‘multirates’), the stiff subsystem with a short step size to achieve stability, and the nonstiff subsystem with a longer step size. These methods with different step sizes are called multirate methods [3], split methods [9], or methods with separate step size [l].
2. The principle of multirate
wo = w(t0)
W,+I = Mw,,
t,,, T; h) .
(2.3)
Here w,+r , w,, are vectors containing all the information necessary to carry out the multistep method. The multirate method can be described with M (2.3): Ynk+i+l
=MCYnk+i,
fnk+i,
T;f),i=o,..-, k- 1 , (2.4)
( bzk+k
= M(bzk,
tnk, kT; S>
with T kT
2 integration step size for the stiff subsystem, 2 integration step size for the nonstiff subsystem, k E N = multirate factor, f:= f(y,0,t)with 0 an approximation of u. A special approximation 0 of the u-values is necessary to evaluate the stiff subsystem f at the points t&+j+l, i = 0, ._., k - 1. At these points no values of IJ are computed by the integration procedure of the nonstiff subsystem. For each n E N the equations (2.4) describe an
methods
If the following n-dimensional lem w’=h(w,t),
with dim@) + dim(u) = n, then the subsystemy’ = f@, u, t) is assumed as the stiff subsystem, Let M be a known linear multistep method applied to (2.1), with the step size T:
initial value prob-
(2.1) 34
W. Gomm /Stability ynk-k
Ynk-k+2
Ynk-k+l
I
analysis of explicit multirate methods
I T
‘.
‘Ynk
I 1
Unk-k
, .’
Ynk+2
Ynk+l
I
T
T
35
I
T “nk+k
“nk
I
I
‘Ynk+k
I kT
kT
Fig. 1. Chronological sequence of the multirate method (points of time of computation).
integration rule for a cycle of the length kT. Therefore the multirate method is regarded as a method with the step size kT (see Fig. 1). If the step size T is given to achieve stability, then the number of the functional evaluation from t,k to tnk+k 1s:
in the normal case: k * funct. eval. off+ in the case of multirate k * funct. eval. off+
k * funct. eval. ofg ,
= min(Pm
of the aims of the stability analysis
The solution of the systems of differential equations which appear by the simulation of dynamical system can be divided into transient and steady-state components. Let w’(t) = h(w, t) + r(t) )
wo = w (to)
method: 1 * funct. eval. ofg .
If the multirate method is used on a serial computer, the evaluation ofg must be spread uniformly over k time steps, in each of them one evaluation off is calculated. Under real-time constraints, this procedure requires very careful programming [4]. Otherwise this problem can be solved on a parallel computer (multiprocessor system with MIMD-structure). In this case one or more processors solve each subsystem [2]. Independent, however, of the programming of the multirate method or of the choice of the computer, it is necessary to investigate the behaviour of error and stability. The order of the local error PMULT depends on the order of the local error pM of the integration method M (2.3) and on the order pi of the approximation 0 ofv: PMULT
3. Specification
PI + 1)
be such a system r s 0 the solution nents, which will For r * 0 and t + state components, function r(t). Let
of differential equations, then for describes the transient compovanish for t + m (stable system). m the solution describes the steadythe system-answer of the forcing
w’=Xw+r,
h E C,
wo=w(to),
r=r(t)
Re(h) < 0, (3.1)
be the test equation for the stability analysis. The general form of a multistep method with 1 steps
f=OliW,_iT,&d,wk-i=O
i=O
(3.2)
applied to (3.1) produces the difference equation:
.
This result is shown in [3] and with some extensions respecting the leading error coefficients in [6]. In this article the stability problem for multirate version of linear multistep methods is described and a way is found to construct the stability polynomial for these methods.
r,_i = r(t,-i)
(3.3)
(the starting values wo, . .. . w[_~ shall be known). The z-transform of the difference equation (3.3) can be found (definition and properties of the z-trans-
W. Gomm / Stability analysis of explicit multirate methods
36
form are described in [7]):
E;f=o piZ’-’
T
‘(‘) =I;&,
so that the steady-state
(cq - XT/$)
R (2) z~-~
G(r):=&e
where W(z) and R(z) are the z-transforms
of the
sequences (wkNand(~&~N. The polynomial in the denominator of (3.4) is called the stability polynomial of the linear multistep method (3.2)
g
(ai - ATPi)Z‘-j =:
p(z) - XTu(z)
.
(3.5)
One of the roots ai, j = 1, .... 1 of the polynomial (3.5) is called the principal roar (here ar), the others are the extraneous roots. It is known that ul is an approximation of eATand that T-’ In ai is an approximation OfX. The method (3.2) is called absolutely stable for AT, if lail < 1 for j = 1, .... 1. \aiJ= 1 is only possible for simple roots aj. The region of absolute stability is the set of all XT, for which the method (3.2) is absolutely stable. The method (3.2) is called rehztively stable for XT, if]a,]>ifzil,j#l. The analysis of the steady-state components of the numerical solution of the differential equation is done only with one of the most important forcing functions: r(f) = eiwf. Using the Laplace-transform the transfer function of the system (3.2) is computed
=-+(s) =-&-&.
W(s)
(3.6)
The factor
A=:
=
is
.
The factor C(iw) = (io - A)-l is called the frequency response, it represents the change of phase and amplitude of the forcing function eiwr, caused by the system. The transfer function of the numerical system is derived from (3.4) (see [5]), it is (with (3.5)):
C(z) =
T@) p(z) - hTu(z) ,=,iwT
so that the steady-state solution is
(3.8)
’
component
of the numerical
Now four aims of a stability analysis can be formulated: (1) to find the region of absolute stability, (2) to investigate the quality of the approximation of ehT by al resp. Xby T-l lnai, (3) to investigate the relative stability, (4) to investigate the quality of approximation of the exact transfer function G(L) by the numerical transfer function C(eiW3. Now it remains to find four analogous points for the stability analysis of multirate methods.
4. Stability analysis for multirate versions of linear multistep methods It is not possible to use the (one dimensional) test evation (3.1) to perform the stability analysis of multirate methods. Therefore at least a two-dimensional test equation is needed:
G(s)
is called the Luplace transfer function of the system. The output function w(t) of (3.6) is w(t)
iwt
component
&
(eiwr - eht)
+&
with the forcing functions with
8
= r(‘)(t), i = 1,2 and
eiWf
(p) let nl, nZ be the eigenvalues of for t+m,
(3.7)
Re(qi) < 0, i = 1,2
LE
with h,-l
(4.2)
W. Comm /Stability
(y) let the first equation
be the stiff equation:
Re(qi) << Re(q2) < 0. The application of the multirate version (2.4) of the general form of a linear multistep method (3.2) to the test system (4.1) leads to the difference equations:
5 j=0
ffjynk+i-i
-
h,T,$
@Ynk+i-j
-
/JT!$
T2
be two different transformation variables for k E N. The z-transform is applied to sequences with the step size kT, the zk-transform to sequences with the step size T. The following theorem shows the transfer from zk- to z-transform: Theorem 4.1 ([7,8]). If the zk-transform of Cjn)na is known:
fljhk+i-j
1 =
Y(zk) =
Pj$A+i-i
c yn zkn ,
with Cj= e2tilk
where the values fin&&j denote the approximations of u,k+j_-j, which are not computed in the second difference equation of (4.3). The difference equations (4.3) define various sequences with different step sizes: the sequence (j,JnE~ with the step size T and its subsequence with the step size kT, the sequence (unk)nCZJJ tink)nEN with the step size kT and a sequence (&k)n~ of approximated values, which has the subsequence (v,&jna. The results of the multirate method (as a method of the step size kT) are two sequences (unk)nEN of the step size kT to approximate tink)ntN> the exact values v(t), v(t) at the points (t,&N. In Section 3 the z-transform method is used to construct the stability polynomial of a linear multistep method applied to the scalar test equation (3.1). To construct the stability polynomial of the multirate method, so-called multirate z-transforms are needed, which are z-transforms with different transformation variables. The problem of these multirate z-transforms is the combination between transforms with different variables. Let
eTs,
then the z-transform of the subsequence
y(z) =; ;ij
(jnk)nEN is
Y(z&‘) (kth root of unity).
The second type of transfer from z- to zk-transform is not unique: it depends on the chosen approximation. The principle of this transfer can be described: let G(Z) = CrZo g&z-” be a known z-transform of the sequence (g&&m. Using suitable approximation methods like inter- or extrapolation polynomials, a piecewise continuous functiong*(t), t Z 0 is constructed, so that
g*(t) continuous on [nkT, (n + 1) kT) for all n E N. Using the Laplace transform method, g*(t) is transformed to L[g,(t)]
= C,(s)
= H(s) G(ekTS) = H(s) G(z)
(see Section 5 for more details). The Laplace transform H(s) is often called interpolator or extrapolator or simply called hold. The Laplace transform G,(S) is transformed to a zk-transform using known tables [7] and the following Theorem 4.2 ([6]). Let the Laplace transform Y(s) = G(z) F(s),
z = ekTs
(-z-transform) be known. If the zk-transform F(zc) of the Luplace transform F(s) can be computed, then
and zk = e Ts
zk =
n=o
for i = 1, .. .. k (4.3)
z = ekTs
31
analysis of explicit multirate methods
= Zl/k
(*
zk-transform)
y(zk)
= G(z)
F(Zk)
>
fk = eTs = zllk
W. Gomm / Srabiliry analysis of explicit multirate methods
38
If the multirate z-transform method is applied to the difference equations (4.3), then an approximation 0 of the u-values has to be found. These approximated values can be computed with extrapolation polynomials from the u- and u’-values at the points t,k-jk, j = 0, 1, ... . In the case of the linear test system (4.1) the following zk-transform can be found, formally expressed: p(zk) = Hh,(zk)
v(z) + H&k)
ycz)
.
(4.4)
In Section 5 the structure and the way of construction of these holds HA,(zk), H,(zk) Will be shown exactly. The following theorem shows the principle of construction of the stability polynomial of the multirate versions of linear multistep methbds.
form to z-transform is only of theoretical importance because the computing effort (c number of complex functional evaluations Y(z&j), j = 0, 1, .... k - 1) becomes larger for larger k. A computing rule with a computing effort independant of k has to be found. 5.1. Construction
of holds
Because of real-time constraints the holds have to fulfil the following conditions: (a) For the construction of holds, only u- and u’values can be used which belong to the past and which are used for the integration of the nonstiff equation. (b) Let t,& be the last used point of time, then the value ukk = g@nk, U&, tnk) cannot be used for the approximation of the values 1, .... k - 1, because the functional evaluation of g needs the whole interval (tnk, tnk+k). Condition (a) allows only the use of extrapolation polynomials for the construction of holds. This can be done by applying the Taylor series, so that the order of error is I- 1: fi,k+j,
Theorem 4.3 ([6]). The stability polynomial of the multirate version (4.3) of a linear multistep method (3.2), applied to the test system (4.1) is: 0 = P(z)
- PTA,(Z)1 [P(Z) - h,kTo(z)l
- /LekTZ u(z) A*,(Z)
(4.5)
ir(nkT+ T) - u(nkT+
with
j
=
T) = O(T~)
(5.1)
for 7 E [0, k7’) and n EN. The set-up for the hold (4.6)
ml O(nkT+ 7) = c
m=O
and
‘5 k B(z):=A,(z)
1
j=O
.
@ktj) P(zkt’)
He(zktj) -
A, T@kf’>
a,(T)
U,k-mk
m2 ,
(4.7)
where 5 = e2nilk is the kth root of unity and p(z), u(z) are from (3.6). k E N is the multirate factor, and the zk-transform of the approximated function 0 is given by (4.4).
5. Implementation of the stability analysis of the multirate versions of Adams-Bashforth methods The application of the Theorem 4.3 to construct a stability polynomial for a concrete multirate method needs some additional considerations: (1) Construction of holds (4.4). (2) Theorem 4.1 about the transfer of zk-trans-
+ T
c m=l
Pm(T)
ULk-mk
(54
leads to the linear system ml
( l&f
.
“1 = gl
(5.3)
[a,(7)(-mIil
+ .EI [Pm(T)(-mP1
forj=l,...,l-1. The approximation i)(t) for t > 0 can be described with the following transformation of variables t=nkT+rer=t-nkT
fornEN:
-J
W.
c c a,?(t n=o m=o
- nkT) u,k-_mk H,(s) = 0 .
m2 +
c m=1
flm(t
.
-
nkT)
t&k-mk
1
[u(t - nkT) - u(t - (n + 1)
(5.4)
Example 5.2. Hold of the order 1 = 1, group B, the set-up (see (5.2)) O(nkT+ 7) = (110(T)u,k + T& (7) t&-k
WI
results in the linear system
with t-co, t>O.
0,
u(t) = { 1,
1 = oo(r) 1
(Let the values u-i, v!_i, j = 1, .. .. max(mr, m2) - 1 be zero because the problem of starting values is not important for the stability behaviour.) The terms am(t - nkT), P,(t - nkT) are polynomials in (t - nkT). Eq. (5.4) is transformed, using Laplace transforms and the linear test equation (4.1)
*
&=;al(T)
P(s) = I&,(s)
V(ekTS) + H,(s) Y(ekTS) .
o(t) = c
n=o
[u,k + (t - nkT) &k-k] . [u(t-nkg--(t
(5.5)
Extrapolation polynomials to the degree 3 are computed in [6], they are divided into two groups: Group A: m2 = 0 in eq. (5.2) Group B: ml = 0 in eq. (5.2). It is H,(s) = 0 for all holds of the group A. Now two examples of holds are computed.
i3(nkT + 7) = 0(,,(r) u,k + @l(T) u,k-k
P(s)= I]
--t--de i
(7)
>
J&Js)
5
=f?+
H,(s) = -St
&=-al(T)
+y)vnk-Funk-k]
. [u(t - nkT) - u(t - (n + 1) k7’)] Application in
of the Laplace transform
kTs 2
)
s2 ekTs
*
Y(ekTs)
p(
_kT+z u( Sz]j sz2
with the solution (ILL = 1 + r/kT, cxl (T) = -r/kT. This solution is set into (5.4):
n=O ((1
ekrs _ 1
kT
so that with z = ekTs:
results in the linear system
1 = (110(T) •t a1
-(n+l)kT)].
Application of the Laplace transform method under consideration of u’(t) = ey(t) + &u(t) results in
.V(ekTS)tE Example 5.1. Hold of the order 1= 1, group A, the set-up (see (5.2))
&(q=j.
The solution is set into (5.4):
161:
i)(t) =
39
so that with z = ekTs:
Ml
m
i)(t) =
Gomm / Stability analysis of explicit multirate methods
method results
z-l __ s2z
’
5.2. Practical execution
of the transfer between zk-
and z-transform
Let g(t) be the time function belonging to the’Laplace transform G(s), then Zk [G(s)] is called the zktransform belonging to the sequence (g(nT)),a. Now the transfer between zk- and z-transform which is used in Theorem 4.3 can be described: The z-transform G(z) of the zk-transform G(zk) has to be found (written: G(zk) o--c G(z)) @k) P@k)
Zk [HAu(s)l -
A, T’J(Zk)
.
_ ‘h”(4 B(z)
W. Gomm /Stability
40
analysis of explicit multirate methods
-=kT
and
zkT
1)2 -m.
(zk-
It follows that The holds in Examples 5 .l and 5.2 and in [6] have a special structure, so that only some fundamental subproblems have to be solved: The calculation of the z-transforms of @k) P(zk)
Zk -
b-j]
for j = 1, 2,....
Tz(~21(4+~2o(~n)) &zkb-21
M
@k)
02r(a,)
(5.6)
’
-a,)-’
+ (1 -a&)
k-‘(1
+ (ah - 1) k-‘(1
-a;)-’
= (1 -a&)-’
,...,
!??=I ,..., 1. (5.7)
The third step requires a lot of considerations will be demonstrated by an example.
which
Example 5.3. Computation
of
&zk[s-21
of the z-transform
m = 1, .... I
7
(see also
~Sj(a,,z)
+----&zk[s-il
the z-transform
Zk [s-‘1 = Tzk(zk - 1)-’ *
1)2 =ZkT
1
1 + (1 -
(1 -a,)2(zk-
&)(Zk
-
1)2+
(1 -
a,)2(zk
- a,)
The following transfers are true [6,7] : zk -•4zk-1
zk -----e----tzk-am
[6])
1
~
Zk ~&&)I
and
Zk-am
-zk[&(~)l Zk-am
can be computed. (2) The decomposition (5.6), (5.7) into partial fractions is remultiplied with the z-transforms mentioned in (1). So the terms (4.6) (4.7) can be computed. (3) Now the polynomial (4.5) can be computed directly. of Example 5.1). The
(z - II2
TZk
@k) =
-1
Tz/c *(zk-a&(zk-
from
of
Example 5.4 (continuation zk-transform of P(s) is
gives [7]:
.
The stability polynomial (4.5) is constructed the separate ‘units’: (1) Using the zk-transforms
1 j=1,2
-u$)-~
and
where 1 is the degree of the denominator (c numbers of steps of the Adams-Bashforth methods). (3) Computation of the z-transform of the terms &zk[s-‘1,
a&)
=: S2(am, z)
(~~~(a~) = (-&)(l
This can be done in 3 steps: (1) Computation of the zk-transform Zk [s-j]. (2) Decomposition of the rational zk-transform into partial fractions
= ii% zk - a,
-
with
.
A, T(-+k)
p(zk) - A,, Tu(zk)
(z - 1)2(z
zk-
1
kT(zk-
1)2 1
v(z)
Z2
1)
1 ’
=: &,(zk)
v(z) .
WithoIl = (1 - ak,)(l -a,)-’ and QZO(U~), a2 1(am) from Example 5.2 the z-transform &
H*u(zk) O---X
[Z(~ll(%)
+ ~21(6n))
Z z-l’
-
Z z-am
k'
(a20(4
-
~ll(4?l))l
Hz
-
dwl
*
Then the stability polynomial of the multirate version of the Adams-Bashforth-2-step method with
41
W. Gomm / Stability analysis of eqolicit multirate methods
the first-order hold of group A is: G22(z)
0 = B(z) [p(z) - x”kTu(z)]
- &W%(z)
=
A (z) with (see also Theorem 4.3)
with
Det(z) = P(z)
A(z)
a!&)
(z -
B(z) = z . *i,
= ?$l bm
- N&(z)1
b(z)
-
Wb~(z)l
- wkT2u(z)A,(Z) >
[z@ll@m)+
Cl(z) =‘(z)
a21(4)
Tmel eiwim -a,
’
C,(z) = kTu(z)
with a,,,, b, from the partial fractions of (5.6). -(a20(~m)-all(am))l
fi
(z-u,3
j#”
where
6 (iw) := [r7, - G(tn)]/rZ(tn)
u(zk) m=l
zk
-am
dzk)
-
In the case of the scalar test equation (see (3.8)), the relative error between the exact transfer function and its numerical approximation is computed by
A, T"czk)
There exists an analogous equation system (5.8)
and
a(&) @k)=
,@k)=&-k,
-l&k+
:= G(eiWkT)[G(io)]-’
- 1.
= G(eiWT) [G(io)]-’ for the transfer
-E,
0.5 .
(5.9)
where
Beside the computation of the stability polynomial the transfer function of the numerical system analogous to (2.8) has to be found. It is a (2 X 2)-system of transfer functions. Using the forcing functions
[(-(ia)]-’
=
i0 - A, -’ --f
iw-
X,
1’
5.3. Numerical results
r(‘)(t) = rc2)(t) = eiWt the transfer matrix can be described by Theorem 5.1 ([6]).
The transfer matrix
1
(5.8)
has the following elements
Gl
=-
1
Det (z)
Gzl (z) = __
1
Det (z)
[P(Z)-
k~“@(Z)l
Cl(Z)
The numerical stability analysis of multirate versions of Adams-Bashforth methods shows a direct connection between the size of the coupling between y and u (=^the elements p, E) and the quality of stability properties of the multirate methods. Therefore a quantitative measurement of the coupling is necessary to describe that phenomena quantitatively. Eq. (4.5) and the computation of the holds show that only the product tie is used in the equations. The size of 1-16 has to be measured in relation to the stiffness ratio, which is given by
9 8B := IRe(rlr)I/IRe(~2)l
PW,(Z)
ekTo(z)
c2 (z) 9
Cl(Z)
7
in the case of the linear equations (4.1) (4.2). The matrices A and a - A, a E R, a > 0, however, must have the same degree of coupling because these matrices have the same stiffness ratio and with the step size T or T/a the same stability polynomial. Therefore the following definition is made:
42
W. Gomm / Stability analysis of explicit multirate methods
Definition 5.2. The degree of coupling GK is measured in relation to the stiffness ratio SR GK := ,t: ,i”,“, %
(sign GK = sign&e)) .
(5.10)
The influence of the size of coupling on the stability properties can be shown in the following way: the multirate method is applied to a set of similar matrices, which are used in (4.3). These similar matrices have the same eigenvalues ql, Q, but different degrees of coupling. The comparative value is the matrix with /_E= 0, that is a degree of coupling GK = 0. In this case the multirate method changes into the application of the linear multistep method to the equation y’ = Xvv = 7,~ with the step size T and to the equation V’= X,LJ= QU with the step size kT. It is well known that the exact solution of the system (4.1) is a linear combination of eqlt and eqzr. But the component e ‘?lr is not a significant part of the solution if the first equation of (4.1) is the stiff equation (see (4.2) 7) and if t > 0. So the second investigation aim in Section 2 is only valid for the principal eigenvalue Q, which is approximated by (kT)-’ In al, where a1 is the principal root of the multirate method. Thus the relative error of the principal eigenvalue can be measured by
The numbers at the curves denote the degrees of coupling which are computed by (5.10). These numerical results show that the multirate versions of the Adams-Bashforth-2-step method have for weakly coupled systems (GK < 5%) nearly the same stability properties like the multirate method applied to the comparative matrix (GK = 0). This is also true for all degrees of coupling, but with step sizes kT < 0.1. To integrate the whole system with the normal version, however, a step size less than 0.02 is necessary to achieve numerical stability. Instead of the error matrix (5.9) the transfer functions of the both componentsy, u and their relative errors are computed to present the numerical results in a clearer way: using the forcing functions r(‘)(t) = rW(t) = eiWr the exact transfer functions can be computed (iw - h,) + P y”(t) =
U(t) =
eiwt =: G,(iw)
eimr ,
eiwf =: G,(iw)
eiwr
(iw - h,)(iw - hY) - PE (iw - X,) + E (io - X,)(io - h,) - ye
and the numerical
transfer functions
(5.11) Some numerical results for the multirate versions of the Adams-Bashforth-2-step method with first order holds and the multirate factor k = 50 are shown in the following figures. One of the hold is out of the group A, the other out of the group B (see the text in the figures). The absolute values of the principal root al (the approximation of e km2) and the absolute value of the absolute greatest extraneous root are shown in the Figs. 2-5(a). The relative errors of the approximation of q2 (computed by (5.11)) are shown in the Figs. 2-5(b). The regions of absolute stability are shown in the Figs. 6,7 (only the upper part of the symmetric regions). In all figures the symbols ‘+’ denote the values, which results from the application of the multirate method to the comparative matrix with GK = PC = 0.
The relative errors can be computed
by
IA&o)1
:= I~Y(eiWk~[Gy(i~)]-l
- 1 I,
IA,(iw)I
:= IG,(e iwkT)[G,(iw)J-’
- 1I .
Table 1 shows the errors of the following example: the multirate versions of the Adams-Bashforth-2-step method with first-order holds, the multirate factor k = 50, and the step size kT = 0.1. The results of this example show that the maximal error of the two components of the normal version differs only a little from the maximal error of the multirate version with the hold of the group A. Both, however, are better than the multirate version with the hold of the group B.
W. Gomm / Stability analysis of explicit multirare methods
43
1.0
0.8
0.6
0.4
0.2
a 0.0
1 0.0
0.1
0.2
0.3
0.4
0.5
0.7
0.6
AB- 2
OROER OF EXTRAP. A 1~ PRINCIPAL ROOT A. MAX. EXTRANEOUROOT Z-PLANE, EIGENtALUESs ( -0.500000E02 , 0.000000E00 ) 1 -O.lOOOOOEOl , o.boooooEoo
0.30 t 0.25
t
AB- 2 ORDER OF EXTRAP. A 1, REL. ERROR OF PRINCIPAL ROOT EIGENhLUES, f -0.500000E02 , 0.OOOOOOEOO ) 1 -0.100000E01
Fig. 2.
,
0.000000E00
I
I
I
0.8
W. Gomm / Stability analysis of explicit multirate methods
44
1.0
0.8
0.6
0.4
0.2
0.0
J 0.0
I
0. I
0.2
0.3
0.4
0.5
0.6
0.7
0.8
AB- 2 ORDER OF EXTRAP. A I, PRINCIPAL ROOT A. MAX. EXTRANEOUSROOT Z-PLANE, EIGENtALUES, ( -0.500000E02 , 0.000000E00 ) ( -O.lOOOOOEOl , O.bOOOOOEOO )
0.30
0.25
0.20
0.15
0.10
0.05
0.00
0.0
0.1
0.2
0.3
0.4
AE- 2 OROER OF EXTRAP. A I, REL. ERROR OF PRINCIPAL ROOT ( -O.lOOOOOEOl EIGENbALUES, I -0.500000E02 , 0.000000E00 )
Fig. 3.
0.6
0.5 ,
0.000000E00
0.7 1
0.8
W. Gomm /Stability
45
analysis of explicit multirate methods
1.0
0.8
0.6
0.4
0.2
0.0
0
)
0.1
0.2
0.3
0.4
0.5
0.6
0.7
AE- 2 ORDER OF EXTRAP. B 1, PRINCIPAL ROOT A. MAX. EXTRANEOUIROOT Z-PLANE, ( -0.100000E01 , 0. boooooEoo EIGENtALUES, ( -0.500000E02 , 0.000000E00 )
)
0.25
AS- 2 ORDER OF EXTRAP. B I, REL. ERROR OF PRINCIPAL ROOT EItENtALUESn I -0.500000E02 , 0.000000E00 I ( -O.lOOOOOEOl
Fig. 4.
,
0.000000f00
)
0.8
W. Gomm /Stability
46
analysis of explicit multirate methods
1.0
0.8 0.4
I 0.6
20
0.4
0.2
a 0.00.0
0.1
0.2
0.3
0.4
0.5
0.6
A&
2 ORDER OF EXTRAP. 6 I, PRINCIPAL ROOT A. MAX. EXTRANEOUIROOT Z-PLANE, EIGEN(rALUES, I -0.500000E02 , 0.000000E00 ) I -0.l00000E0l , o.booooofoo
0.7
0.8
)
0.30
0.25
0.20
0.15
0.10
0.05
0.00
0.0
0.1
0.2
0.3
0.4
0.5
AS- 2 OROER OF EXTRAP. S 1, REL. ERROR OF PRINCIPAL ROOT ( -O.lOOOOOEOl EIGENtALUES, I -0.500000E02 , 0.000000E00 )
F3g. 5.
0.6 ,
0.000000E00
0.7 )
0.8
W. Gomm / Stability analysis of explicit multirate methods
41 1.0
a
AB- 2,
ORDER
OF
EXTRAP.
A
I,
REGION
OF STAB.,
STIFFNESS
RATIO
50
1.c
AEI- 2, ORDER
OF
EXTRAP.
A
1,
REGION
OF STAB.,
Fig. 6.
STIFFNESS
RATIO
50
W. Gomm / Stability analysis of explicit multirate methods
48
1.0
-1 AB-
:o
2, ORDER
OF EXTRAP.
Bl,
REGION
OF STAB.,
STIFFNESS
RATIO
50
1.0
b
AB-
2, ORDER
OF EXTRAP.
01,
REGION
OF STAB.,
Fig. I.
STIFFNESS
RATIO
50
W. Gomm / Stability analysis of explicit multirate methods
Table 1 Adams-Bashforth-2-step GK
+I%
ilO%
?50%
w
method,
k = 50, kT=
Normal
49
0.1
version
1. order hold B
1. order hold A
IA,(iw) I
I A,(iw)
0.0 0.0
0.33E-1 0.33E-1
0.80E-3 0.82E-3
0.85E-3 0.83E-3
0.47E-3 0.47E-3
0.35E-1 0.35E-1
0.73E-3 0.1 lE-2
0.29E-2 0.29E-2
0.35E-2 0.34E-2
0.29E-2 0.29E-2
0.37E-1 0.37E-1
0.25E-2 0.35E-2
0.25E-5 0.24E-5
0.35E-1 0.35E-1
0.21E-1 0.21E-1
0.35E-1 0.35E-1
0.43E-1 0.43E-1
0.35E-1 0.36E-1
0.33E-5 0.32E-5
O.lOEO O.lOEO
0.41E-1 0.40E-1
O.lOEO O.lOEO
0.50E-1 0.50E-1
O.lOEO O.lOEO 0.77E-2 0.84E-2
IA,(b) I
Ih&w) I
IA&w)
0.0
0.0 0.0
0.0 0.0
0.0 0.0
0.5
0.77E-7 0.75E-7
0.47E-3 0.47E-3
1.0
0.43E-6 0.43E-6
3.0
5.0
I
I A,(k)
I
-
0.0
0.0 0.0
0.0 0.0
0.0 0.0
0.0 0.0
0.60E-1 0.52E-1
0.5
O.l3E-6 O.l3E-6
0.47E-3 0.47E-3
O.l5E-2 O.l3E-2
0.49E-3 0.50E-3
0.62E-1 0.5 7E-1
0.73E-2 0.84E-2
1.0
0.81E-6 0.78E-6
0.29E-2 0.29E-2
0.67E-2 0.62E-2
0.30E-2 0.30E-2
0.69E-1 0.67E-1
0.5 7E-2 0.1 OE-1
3.0
0.72E-5 0.67E-5
0.35E-1 0.35E-1
0.5 3E-1 0.50E-1
0.35E-1 0.36E-1
O.llEO 0.89E-1
0.27E-1 0.43E-1
5.0
0.14E-4 O.l2E-4
O.lOEO O.lOEO
0.12EO O.llEO
O.lOEO O.lOEO
0.14EO 0.14EO
0.94E-1 O.llEO
0.0
0.0 0.0
0.0 0.0
0.0 0.0
0.0 0.0
0.94E-1 0.36E-1
0.37E-1 0.44E-1
0.5
O.l6E-6 O.l5E-6
0.47E-3 0.466-3
0.24E-2 O.l3E-2
0.94E-3 0.97E-3
0.95E-1 0.50E-1
0.36E-1 0.43E-1
1.0
O.lOE-5 0.97E-6
0.29E-2 0.29E-2
0.94E-2 0.70E-2
0.38E-2 0.41E-2
O.lOEO 0.73E-1
0.32E-1 0.43E-1
3.0
0.1 lE-4 O.lOE-4
0.36E-1 0.35E-1
0.80E-1 0.74E-1
0.35E-1 0.39E-1
0.15EO OS7EO
O.llE-1 0.76E-1
5.0
0.26E3-4 0.23E-4
O.lOEO O.lOEO
0.20EO 0.18EO
0.97E-1 O.llEO
0.23EO 0.25EO
0.63E-1 0.15EO
6. Discussion The necessity for a two-dimensional test equation is the essential point of difference between this analysis and the well-known stability analysis of the ‘normal’ multistep methods. We have shown a general construction rule for the stability polynomial of the multirate versions of linear multistep methods, using the
I
multirate z-transform method. This rule has been applied to the stability analysis of multirate version of the Adams-Bashforth methods and has been demonstrated in detail on the Adams-Bashforth-2-step method. The stability properties of these multirate versions differ for weakly coupled systems slightly from the stability properties of the methods applied to the comparative matrices.
50
W. Gomm /Stability
analysis of explicit multirate methods
References [l]
H. Blum, Numerical integration of large scale system, using separate step size, in: A.W. Bennett and R. Vichnevetsky, Eds., Numerical Methods for Differential Equations and Simulation (North-Holland, Amsterdam, 1978) 19-23. [2] W. Bub, Private communication. [ 3) C.W. Gear, Multirate methods for ordinary differential equations, Department of Computer Sciences, University of Illinois, Report UIUCDCS-F-74-880. [4] C.W. Gear, Simulation: conflicts between realtime and software in: J.R. Rice, Ed., Mathematical Software III (Academic Press, New York, 1977) 121-138.
[ 5I E.G. Gilbert, Dynamic-error
analysis of digital and combined analog-digital computer systems. Simulation 6 (1966)241-257. [6] W. Gomm, Stabilit&suntersuchung van expliziten Multirate-Methoden zur numerischen Liisung von Anfangswertproblemen bei gewljhnlichen Differentialgleichungen, Dissertation, TU Braunschweig (1980). of the z-Transform (7 E.I. Jury, Theory and Application Method (Wiley, New York, 1964). ]8 I M.A. Neighbors and CA. Halijak, Effective computation of exact multirate z-transforms, in: Vogt and Mickle, Eds., Modelling and Simulation Vol. 8, Past 1, Proc. 8th Annual Pittsburg Conference (1977) 539-542. methods for simultaneous [9 J.R. Rice, Split Runge-Kutta equations, J. Res. Nat. Bureau of Stand. 64B (1960) 151-170.