CC BY
36
MSC 65L80, 65L05, 65L06, 65L20
DOI: 10.14529/ mmp180406
CONVERGENCE ANALYSIS OF LINEAR MULTISTEP METHODS FOR A CLASS OF DELAY DIFFERENTIAL-ALGEBRAIC EQUATIONS
Vu Hoang Linh1, Nguyen Duy Truong2, M.V. Bulatov3
1 Faculty of Mathematics, Mechanics and Informatics, Vietnam National University,
Thanh Xuan, Hanoi, Vietnam
2
3
Russian Academy of Sciences, Irkutsk, Russia
E-mails: [email protected], [email protected], [email protected]
Delay differential-algebraic equations (DDAEs) can be used for modelling real-life phenomena that involve simultaneously time-delay effect and constraints. It is also known that solving delay DAEs is more complicated than solving non-delay ones because approximation of solutions in the past time is usually needed and discontinuity in higher derivatives of the solutions is typical. Recently, we have proposed and investigated linear multistep (LM) methods for strangeness-free DAEs (without delay). In this paper, we extend the use of LM methods to a class of structured strangeness-free DAEs with constant delay. For the approximation of solutions at delayed time we use polynomial interpolation. Convergence analysis for LM methods is presented. It is shown that, similarly to the non-delay case, the strict stability of the second characteristic polynomial associated with the methods is not required for the convergence if we discretize an appropriately reformulated DDAE instead of the original one. Numerical experiments are also given for illustration.
Keywords: delay differential-algebraic equation; strangeness-free; linear multistep method; stability; convergence.
Dedicated, to Professor Viktor Chistyakov on the occasion of his 70th birthday.
Introduction
Delay differential-algebraic equations (DDAEs) arise as mathematical models of real-life processes in which time lag and constraints simultaneously appear. Areas of applications include electrical circuit design, real-time simulation of mechanical systems, chemical engineering, power systems, control and optimal control, etc, see [1-4]. While the theory and numerical analysis of delay (ordinary) differential equations (DDEs) as well as those of DAEs (without delay) have been fairly well established, see [2,5,6] and [7-9], respectively, the same cannot be said in the case of delay DAEs which analytical and numerical solutions have not been completely understood yet. Even the solvability of general linear delay DAEs has been investigated only very recently in [3,10]. Very few papers have been devoted to the convergence analysis of numerical methods for delay DAEs and most of them are restricted to consideration of delay DAEs in semi-explicit form and implicit numerical schemes. In [11], Ascher and Petzold investigated BDF and collocation Runge-Kutta methods for semi-explicit DDAEs of retarded and neutral type with single delay. Later, Hauber extended the use of collocation methods to retarded DDAEs of index two with state-dependent delay [12]. Liu and Xiao studied the convergence of linear multistep and one-leg methods for semi-explicit index 2 DDAEs with variable delay [13].
Difficulties that arise in solving DDAEs were discussed in [1,2,4]. It was pointed out that in general DDAEs are neither DAEs nor DDEs. However, under certain appropriate conditions, a DDAE may be reduced to a DDE of retarded or neutral type [1].
Recently, in [3,10] Ha et.al. investigated general linear variable coefficient DDAEs with constant delay
E(t)x'(t) = A(t)x(t) + B(t)x(t - t) + y(t),
(1)
which may arise, for example, as the result of linearizing general DDAE F(t, x'(t), x(t), x(t — t)) = 0 around a reference trajectory. They proposed an algorithm for regularization of DDAE (1), i.e. a procedure that reduces the original system to a regular strangeness-free DDAE of form
Ei(t) 0
x' (t)
Mt)
Mt).
x(t) +
Bi(t)
Mt).
x(t - t) + î(t),
(2)
where
-EJl(t)
Mt)
methods can
is pointwise invertible. Then, it was also shown that s-stage collocation
i)e implemented for (2) and they are convergent of order at least s. In this paper, we consider a more general class of nonlinear structured DDAEs of form
f ht,x(t),x(t — t ),E (t)x'(t)) = 0, ght, x(t),x(t — t)) =0
(3)
for all t G I = [0, T], t > 0 is a constant delay, which clearly includes DDAEs (2) as a special case. Here we assume that E G C 1(I, R™1™) and fonctions f (t,u,v,w) : I x Rm x Rm x Rmi ^ Rmi, g(t, u,v) : I x Rm x Rm ^ R™2 ,m1 + m2 = m are sufficiently smooth functions with bounded partial derivatives. Given initial condition
x(t) = ф^) for t e [-t, 0],
(4)
we suppose that initial value problem (IVP) (3), (4) has unique solution x(t). Here, x is said to be a solution if the followings hold: (i) It is continuous and piecewise continuously differentiate on I; (ii) It satisfies DDAE (3) for t G I pointwise except for a finite number of discontinuity points as well as initial condition (4). In this paper, we assume that Jacobian
fw E (t)
9u
is nonsingular
(5)
x(t)
x
h)
function 0 must be consistent, that is g(O,0(O),0(—t)) = 0 and 0 £ C([—t, 0], Rm) at least.
The main aim of this paper is to extend the use of linear multistep methods that are proposed for non-delay DAEs in [14,15] to DDAEs (3). Similarly to the approach used in [15-17], instead of direct discretization, numerical schemes are applied to reformulated form
f(t,x(t),x(t - t), (Ex)'(t) - E'(t)x(t)) = 0, g(t, x(t),x(t - t)) = 0.
(6)
In the implementation, approximation of numerical solutions at retarded time may be needed and it is done by using interpolation or continuous extension. The linear multistep methods proposed in this paper provide an alternative approach in addition to the existing BDF and collocation methods. Moreover, discretizations based on explicit methods require less computational cost than implicit ones in the case of large-sized and nonstiff problems.
The organization of the paper is as follows. In Section 1, we present some preliminary results including the method of steps and the analysis of DDAEs (3) by using transformation and reduction. Construction and convergence analysis of linear multistep (LM) methods combined with polynomial interpolation are given in Section 2. In Section 3, some numerical experiments are carried out to illustrate the theoretical results in previous sections. We close the paper by conclusions in Section 4.
1. Preliminary 1.1. Method of Steps
For DDEs with constant delay, the method of steps is a standard tool to investigate analytical as well as numerical solutions [5]. Analogously, this method can be extended to the analysis of DDAEs. Here I VP (6), (4) is replaced by a sequence of the IVPs on the time intervals [It, (l + 1)t] for nonnegative integer / provided that x is known on interval [(/ — 1)t, It]. Therefore, we obtain a sequence of "local" IVPs for non-delay strangeness-free DAEs of the form
f (t,xi+i(t),xi (t — t), (Exi+i)'(t) — E'(t)xl+i(t)) = 0,
U u\ u W n t G [lT, (l + 1)t] ^ g[t,xi+i(t),xi (t — t )) =0,
together with initial conditions
xi+i(/T ) = xi(/T), l = 0,1,... (8)
For l = 0, we define x0(t) = $(t — t), 0 < t < t. Assuming that all the initial conditions are consistent, i.e. g(j,T,xl(lT),xl((l — 1)t)) = 0, the solvability of IVPs (7), (8) for all l = 0,1,..., implies the solvability of original I VP (6), (4). Then, we set
x(t) = xl+i(t) if t e [It, (l + 1)t], l = 0,1,....
Clearly, the "global" solution x(t) is continuous and piecewise continuously differentiable. At the connecting points It, discontinuity in the first or higher derivatives of x is typical. For both DDEs and DDAEs with a single constant delay t, the discontinuity happens at points It, / = 0,1,..., see [1]. Moreover, the existence, uniqueness and smoothness of solutions depend on given initial function $(t). Throughout this paper, we assume that initial function $(t) is consistent and sufficiently smooth such that the unique solution of IVP (3)-(4) is continuous on [0, T] and sufficiently smooth on each subinterval [It, (l + 1)t].
1.2. Reformulation and Conditioning
Conditioning analysis of semi-explicit DAEs of index less or equal to two was considered first in [18]. Then, the authours extended the analysis to semi-explicit DDAEs of retarded and neutral type with single delay in [11]. They showed that a semi-explicit DDAE of index
1 is well conditioned when the essential-underlying-delay ODE (EUDODE) associated with the DDAE is well conditioned. For implicit DDAEs like (3) or even more general ones, the difficulty relies in the definition of EUDODE and the classification of the problem because differential and algebraic variables are not separated as in the case of semi-explicit DDAEs.
Returning to the structured strangeness-free DDAEs (3), by the same approach as in [15, 17], the special structure of (3) is exploited to transform DDAEs (3) into reformulated form (6) which are equivalent to semi-explicit index-1 DDAEs. By condition (5), there exists a pointwise invertible matrix function Q(t) = [Q(1) Q^], where Q(1 £ C 1(l, R™), Q(2) £ C 1(I, R™), such that EQ = [I 0]. Thus, we can assert that
matrices fw, guQ(2\ E are nonsingular. Using the change of variables x = Qy =
[gu]
Q(1'ly1 + Q(2)y2, we have that (Ex)'(t) = (EQy)'(t) = y[(t). Therefore, system (6) becomes
f ht, Q(t)y(t),Q(t — t)y(t — t),y[(t) — E'(t)Q(t)y(t)) = 0,
g(t,Q(t)y(t),Q(t — t )y(t — t ))=0.
By the Implicit Function Theorem, there exists a function f such that from the first equation of (9), we have
y[(t) — E' (t)Q(t)y(t) = f{t,y1(t),y2(t),y1(t — t ),y2(t — t )). (10)
Then, the system (9) is rewritten as
y\(t) = f{t, y1 (t),y2(t),y1 (t — t),y2(t — t)), 0 = g{t,y1 (t), y2(t),y1(t — t),y2(t — t)),
(11)
where f(t,y1(t),y2(t),y1(t - r),y2(t - r)) = f(t,y1(t),y2(t),y1(t - r),y2(t - r)) + E'(t)Q(t)y(t), g(tJy1(t),y2(t),y1(t-r),y2(t-r)) = g(t,Q(t)y(t),Q(t-r)y(t-r)).Itiseasy to check that = guQ(2 which is nonsingular. Hence, system (11) is a semi-explicit
index-1 DDAE of retarded or neutral type that is analyzed in [11]. All the discussions on the conditioning of the IVPs for (11) that was studied in [11] can be applied. Again by the Implicit Function Theorem, there exists a function g such that the algebraic variable of (11) is expressed in the form
y2(t) = g(t,yi(t),yi(t - T),y2(t - T)). (12)
We consider two cases.
a) If y2(t - t) does not appear in the second equation (11), then the equation (12) becomes
y2(t)= g(t,yi(t),yi(t - t)). (13)
t - T
y2(t - t) = g(t - T,yi(t - t),yi(t - 2t)). (14)
Inserting (13), (14) into the first equation (11), we obtain a DODE of retarded type
y'i (t) = 7( t, yi (t) ,g(t,yi (t) ,yi(t - t )),yi (t - t ),g(t - t, yi (t - t ), yi (t - 2t ))), (15)
which can be also written by renotation as
y[(t) = /(t,Vi(t),Vi(t - T),y1(t - 2t)).
(16)
The DDAE (3) is well conditioned (in the sense that its solutions are not too sensitive to small changes in the initial function) if EUDODE (16) is well conditioned and transformation Q is well conditioned.
b) In the general case, we have to propagate the recursion in (12) back from ¿to t — It, where —t < t — It < 0, i.e. t E [It, (l + 1)t] is assumed. For simplicity, we suppose that linearization is applied, equation (12) is written in the form
y2 (t) = g( t,yi(t),yi (t — t ))+ R(t)y2(t — t ), where R(t) = — {guQ(2)(t))-1 gvQ(2)(t — t). This gives
(17)
—1
V2 (t)={H R(t - jt )] V2(t - It ) + j=0
l-1 i-1
+ Z [ n R(t - JT^ g{t - IT, V1 (t - IT),yi(t - IT - T)).
i=0 j=0
(18)
A similar formula holds
l-1
V2(t-T )= n R(t-JT ) V2(t - lT ) +
j=1
l 1 [ i 1
(19)
+ Z [ n R(t - JT )] 9{ t - iT,Vi(t - IT ),Vi(t - IT - T ) )
i=1 j=0
l
y'l (t) = F(t,yi (t) ,yi(t — t ),yi(t — 2t ),..., yi(t — It ),y2(t — It )),
(20)
which is actually DDE of a neutral type. The well-conditioning of DDAEs (11) depends not only on this DODE but also on factor R. If supi>0 ||R(t)|| < 1, Q is well-conditioned, and DODE (20) is well conditioned, then DDAE (11) is well-conditioned, too. The above analysis is demonstrated in the following example.
Example 1. We consider strangeness-free DDAE with constant delay t > 0 of form
1 -ut 0 0
x' (t)
X u(1 - Xt) -1 (1+ut)
+
0 a
b c - bu(t - t)
x(t - T) -
x(t) +
aeMt-r)
(b + c)ex(t-T \
(21)
on interval [0, T] with real parameters a, b, c, u, X. System (21) possesses analytical solution
" ext (1 + ut)
x
.,\t
provided that the initial function is given equal to the exact solution for all —t < t < 0.
If we take Q(t)
becomes
1 ut 0 1
using the change of variables x(t) = Q(t)y(t) system (21)
10 00
y ( t) =
A 0
11
y(t) +
0 a b c
y(t - t) -
aex(t-Tj (b + c)ex(t-Tj
(22)
where y(t) = [y1(t),y2(t)]T. System (22) is rewritten as
y[(t) = Xy1(t) + ay2(t — t ) — aex(t-T
y2(t) = —cy2 (t — t ) + y1(t) — by1(t — t ) + (b + c)ex(t-T \
If c = 0, we obtain algebraic constraint
y2(t) = y1(t) — by1(t — t ) + (b + c)ex(t-T and DODE of retarded type with two lags
y[(t) = Xy1(t) + ay1(t — t ) — aby1 (t — 2t ) — a(b + c — l)ex(t-T \
(23)
(24)
If DODE (24) is well conditioned and uT is of moderate size, then DDAE (21) is well conditioned, too. This property depends on the choice of parameters A, a and b. If c = 0, the second equation (23) leads to
y2(t) = cly2(t - It ) + H( t,yi(t),yi(t - T),..., yi(t - It )):
(25)
where H is an appropriately defined function and —t < t — It < 0. Then, we also have
y2(t — T ) = cl-1y2(t — It ) + Hht — T,y1(t — t ),y1(t — 2t ),..., y1(t — It )). (26) Substituting (26) into the first equation (23), we obtain DODE with l lags
y\(t) = Ay1 (t) + H ht,y1(t — t ),y1(t — 2t ),..., y1(t — It )) + acl-1y2 (t — It ), (27)
which is actually neutral DODE. Therefore, we say that DDAE (21) is of neutral type. If
c
satisfying \c\ < l or \c\ > l, but of moderate size, then problem (21) is well conditioned. Fig. 1 shows actual errors, the exact solution (Exact. Sol) and a numerical solutions (Num. Sol) of well-conditioned problem (21) with time-delay t = l and parameter set A =
— l, 5,u = 10, a = 0, 5,b = l,c = 0, 8. The results of ill-conditioned problem (21) with time-delay t = l and parameter set A = —l, 5,u = l0,a = —2,b = l, 5,c = l, 2 are plotted in Fig. 2. For both problems, time interval [0, 20] is used.
For retarded DDAE (21) with t = l, the results of well-conditioned problem with A =
— l,u = l,a = — l, b = l, 5,c = 0 on time interval t £ [0, lO] are plotted in Fig. 3. It should be noted that a small delay might make a well-conditioned problem be ill-conditioned. This is demonstrated in Fig. 4 for problem (21) with A = —l,u = l, a = —l,b = l, 5,c = 0 and t = 0, 2.
For the above illustrations, the numerical solutions are computed by 2-step half-explicit Adams Bashforih (HEAB2) method, which is presented in the next section.
Fig. 1. Numerical results for well-conditioned problem (21)
Fig. 2. Numerical results for ill-conditioned problem (21) 2. Linear Multistep Methods
In this section, we construct and analyze half-explicit linear multistep (HELM) methods for reformulated DDAEs (6). Consider linear multistep (LM) method which coefficients ai, Pi, (i = 0,1,..., k) satisfy a0 = 0 and [s = 0 where [s denotes the first non-zero coefficient among [^s. If s = 0, then the method is said implicit. Otherwise, it is explicit. LM scheme applied to I VP for ODEs
y'(t) = X(t,y), y(to) = Уo,
is proposed as follows
^ ^ aiyn-i ^^ ^ [ix(tn—il yn-i) ■ i=0 i=s
(28)
(29)
Fig. 3. Numerical results for well-conditioned problem (21)
Fig. 4. Numerical results for ill-conditioned problem (21)
Construction and numerical analysis (order, zero-stability, convergence) of linear multistep methods for ODEs are given in details in numerous books, e.g.. see [7,19].
We assume that the unique solution of I VP (3), (4) is continuous on [0, T ], p times continuously differentiate and (p + l)-th derivative is bounded on each subinterval [It, (l + 1)t]. This means that at points It higher derivatives of the solution may not exist. Assuming w.l.g that the length of integration interval T is multiple of t, we introduce a mesh n = {0 = t0 < t1 < ... < tN = T} which include all "discontinuity points" It. The mesh is not necessarily uniform, but just for the sake of simplicity we assume that a uniform mesh with stepsize h is used. It is possible to extend the construction of LM schemes and their analysis to the case of variable stepsize. Let us denote by xn, Wn the approximations of exact solution x(tn) and derivative W(tn) := (Ex)'(tn). Suppose that approximate solution obtained at the mesh points of n is denoted by xh. For implementation, we will solve numerically IVP (3), (4) on successive intervals [It, (l + 1)t] for any nonnegative
integer l. At time step t = tn, (It < tn-k < tn < (l + 1)t), we assume that previous approximations xn-1,xn-2,... ,xn-k, and Wn-2, Wn-3,..., Wn-k are given. Moreover, we need to approximate the values of x at retarded times tn — t and tn-s — t which may not be mesh points (for example, if the mesh is not uniform). If l = 0 then those values are given by the initial function. Otherwise, they are approximated by local interpolant . through values of xh at mesh points.
From equation (29), we have
1 k k Pi
Wn-s = T-r y^aiEn-ixn-i — y] Wn-i, (30)
hpi 1=0 i=S+i P1
where En-i = E(tn-i). A LM scheme applied to reformulated DDAEs (6) is proposed as follows
f {^k'n-s, xn-s, px (tn— S T)) hp ^ y aiEn-ixn-i ^y ^ p ^Wn-i E n-s xn-^ 0,
kk
hp ^aiEn-ixn-i — ^ p ■■■n-% - US-US J , .
hPs i=0 i=s+1 Ps y61)
g(tn,xn,pxh(tn — T)) = 0,
here E'n-s = E'(tn-s). System (31) is equivalent to
kk
f ^tn-sj xn-s, px (tn—s T) , ^ ^ an En-ixn-i ^ ^ ps E n-sxn-s^ 0j
n s n s n s- n i n i- n i- n s n s
1=0 hps ^^s > (32)
g(tn,xn,pxh(tn — T)) = 0.
We obtain a nonlinear equation for unknown variable
xn
follows
Hn(x,n,x,n-1,... ,xn-k, Wn-2, Wn-3,..., Wn-k, pxh(tn — T ),.xh(tn-s — t ), h) = 0, (33) where Jacobian matrix of Hn with respect to
xn
d
dx,
Hn (xnj xn- ... j xn-k j Wn-2, Wn-3, . . . , Wn-k, pxh(tn — T), pxh(tn-1 — T ),h)
n
(k k ¡3- \
tn-sj xn-si .x (tn— s T) , y ] h3 En-ixn-i 33~^^n-i E n-sxn-s J En
i=0 i=s+1 J
gu (tn,xn,pxh (tn — T)) ,
(34)
h
xn h x
xn
be approximated by Newton's iterative method. After that, approximation Wn-s that will be used for the next step is computed by (30).
We note that the computational procedure should be implemented carefully. First, we have to evaluate the starting values on each subinterval [It, (l + 1)t]. At time t = tn, approximations <.pxh(tn-s — t) and .xh(tn — t) are interpolated by using approximate values at mesh points close to tn-s — t and tn — t, respectively, provided that these node points must belong to interval [(l — 1)t, It]. Therefore, formulas of interpolant . depend on the
concrete location of tn—s — r and tn — r, but they are assumed to have the same order. If the underlying LM method is explicit (s > 1), then scheme (31) for DDAEs (31) is called half-explicit. Otherwise, it is said an implicit LM scheme for DDAEs (31). In the following, we analyze the convergence of LM scheme (31).
Theorem 1. Consider well-conditioned IVP (3), (4) and k-step LM method, which is zero-stable and of order km > 2 for ODEs. We assume that interpolant ф is accurate of order O(hka), ka > 2. Then, k-step LM scheme (31) is convergent of order q = min(p,km,ka), i.e.,
max \\x(tn) — xn\\ = O(hq),
0<n<N
provided that the starting values are accurate to O(hq)
Proof. We use the method of steps and prove the convergence of LM method (31) on each subinterval [lr, (l + 1)т]. First, we consider interval [0, r] and time step tn (0 < tk < tn < r). We suppose that startin g values x0,x1,... ,xk—1, and W0 ,W1,... ,Wk—s—1 are accurate to O(hq). The retarded values are given by ф^п—,3 — т) = — r), and
ф(^ — т) = Ф^п — r). Therefore scheme (31) becomes
1 k k в' f ( tn— S , xn-s , ф(tn—s т^ , 7 Q ^ ^ aiEn-ixn-i ^ ^ Q E n-sxn-sj 0 , , 4
hes i=0 i=s+1 es
g{tnjXnjф(tn — r)) = 0.
This is exactly LM method applied to DAE (7) that is analyzed in [15]. By [15, Theorem 4], we conclude that scheme (31) is convergent of order q on [0, r], i.e., \\x(tn) — xn\\ = O(hq)
0 < tn < т.
Next, we consider interval [lr, (l + 1)r] for l > 1, and assume that approximates xi = x(ti) + O(hq), (ti < lr) are given. We assume that start ing values xmi+1,xmi+2,..., xmi+k—1 and Wmi+1, Wmi+2,..., Wmi+k—s— 1 are already given and accurate of order q, where tmi = lr. At time t = tn (ml + k < n < mi+1), from the assumption of interpolant ф, we have
фXh(tn — r ) = x(tn — r) + O(hq), фXh(tn—s — r ) = x(tn—s — r) + O(hq). (36)
Substituting (36) into system (31) yields
1 k k ei
f ( tn—s , xn—s , x(tn—s т) + O(h ) , , Q ^У ^ aiEn—ixn—i ^У ^ n ^Wn—i E n-SXn-Sj 0, /q>7\
hes i=0 i=s+1 es
g(tn,xn,x(tn — r) + O(hq)) =0.
It follows that
1 k k в i hPsH tn—s, xn—s, xl (tn—s т^ , 7 Q ^ ^ aiEn—ixn—i ^У ^ Q E n—sxn—sj ^nj iOQ \
hes i=0 i=s+1 es
g (tni xni xl (tn т^ 9nj
where 5n = O(hq+1), 9n = O(hq). This is exactly perturbed LM scheme applied to DAE (7), whose analysis is already studied in [15]. Let us denote by {xn} the solution of following (unperturbed) LM scheme applied to DAE (7)
k k
hPs i=0 i=s+1 ' (39)
g(tn,Xn ,Xl(tn - T)) = 0,
hfisf^tn-s,xn-s,xl (tn—S T), hps ^ y aiEn-iXn-i y ^ p ^^n-i E n—s Xn—s ^
where the starting values are taken by exact values x(tmi+1),x(tmi+2),... ,x(tmi+k—1) and W(tmi+i),W(tmi +2),..., W(tmi+k-s-i)■ By [15, Theorem 2], we conclude that \\Xn - xn\\<C max \\x(ti) - Xi|| + D max \\Wt) - Wi\\ +
mi<i<mi +k—1 mi<i<mi+k—s—1 . .
+ K max \\5i/h\\ + L max \\0i\ = O(hq) ^ '
mi+k<i<mi+i mi+k<i<mi+i
for all ml + k < n < ml+1, where C, D, K, L are constants independent of h. From [15, Theorem 4], it follows that
\\x(tn) - Xn\\ = O(hq) for all ml + k < n < ml+1.
By combining the above estimates, we have
\\x(tn) - xn\\ < \\x(tn) - Xn\\ + \\Xn - xn\\ = O(hq) (41)
for all ml + k < n < ml+1. Therefore, we conclude that the global error of x on interval [It, (l + 1)t] is O(hq). The analysis is repeated on successive interval [(l + 1)t, (l + 2)t] in the same way. Finally, by induction, we obtain the global convergence of scheme (31) on interval [0, T].
□
Remark 1. The implementation of LM methods is more complicated than in the non-delay case because we must calculate starting values on each subinterval [It, (l + 1)t]. This extra cost can be reduced when the solution of I VP (3), (4) is sufficiently smooth globally on [0, T]. Then, LM scheme (31) is applied right from n = 1.
Furthermore, if we use a uniform mesh n* with stepsize h = M for some integer M > 1, then retarded values x(tn - t) and x(tn—s - t) are taken by xn—M and xn—M—s, respectively. Therefore, LM scheme (31) is rewritten as
1 k k s
fitn—s,xn — s ,xn — M —s,^Q y ^ aiEn—ixn—i y ^ Q ^Vn — i E n—sxn—sj 0, //(O^
hps i=0 i=s+1 Ps l42)
g{tn,xn,xn—M) = 0.
We do not need to call any explicit interpolation at each step of the computation. It is easy to verify that in this case LM method is convergent of order q = min(p, km).
Remark 2. For more general DDAEs of form
f (t,x(t),x(t - T),x'(t)) =0, g(t,x(t),x(t - t)) =0,
direct discretization by LM methods can be realized analogously to the non-delay case that is proposed and analyzed in [14]. However, for the stability of this direct discretization, an extra condition is required, namely the second characteristic polynomial associated with underlying LM method must be strictly stable, i.e. all their roots must lie inside the unit complex disk. For example, among the popular classes of LM methods, Adams-Moulton methods do not fulfill this condition.
3. Numerical Experiments
We consider again Example 1 with time-delay r =1. We first check the convergence of the half-explicit two-step Adams Bashforih (HEAB2) method that has already been used in Section 1. System (21) is reformulated as follows
[xiit] — utx2(t))' = Xx\(t) — Autx2(t) + ax2(t — 1) — aex(t-1) xi(t) — (1+ ut)x2(t) = bxi(t — 1) + (c + bu — but)x2(t — 1) — (b + c)ex[t-1). ^ '
Discretizing (44) by HEAB2 method on a uniform mesh and using 4-node forward Newton interpolation, numerical results of two well-conditioned problems are displayed in Tables 1 and 2. Here, we compute actual errors ei(h) = max \xi(tn) — xi>n\, i = 1, 2 for various
step-sizes. The numerical convergence order is also estimated by rate = log2 ^)) • The computations are implemented in Matlab. The second-order convergence of HEAB2 method is confirmed in the case of well-conditioned IVPs. However, Figs. 2 and 4 in Section 1 show that it is not the case with ill-conditioned problems, for which the errors
t
Table 1
Numerical results for (21) on interval [0, 20] with A = —1, 5,u = 10, a = 0, 5, b = 1, c = 0, 8
h = 0,03 HEAB2 method
Error in x1 x1 Error in x2 Rate in x2
h 6,9380e-03 _ 3,4484e-04 _
h/2 l,7201e-03 2,0120 8,5222e-05 2,0166
h/4 4,2736e-04 2,0090 2,1173e-05 2,0090
h/8 l,0650e-04 2,0047 5,2760e-06 2,0047
h/16 2,6580e-05 2,0024 l,3168e-06 2,0024
h 32 6,6394e-06 2,0012 3,2893e-07 2,0012
Table 2
Numerical results for (21) on interval [0, 5] with A = -2, ш = 1, a = -2, b = -1, 5,c = 1, 5
HEAB2 method
h = 0,03 Error in x1 x1 Error in x2 Rate in x2
h 9,7882e-04 _ 5,7463e-04 _
h/2 2,4387e-04 2,0050 l,4062e-04 2,0308
h/4 6,0642e-05 2,0077 3,4811e-05 2,0142
h/8 l,5107e-05 2,0051 8,6617e-06 2,0068
h/16 3,7692e-06 2,0029 2,1604e-06 2,0033
h 32 9,4129e-07 2,0015 5,3949e-07 2,0016
We perform two other experiments with first well-conditioned problem (21). The numerical results are obtained by 3-step HELM (HELM3) method and they are compared with those by the standard HELM3 method (i.e. the direct discretization) in Table 3. Here, the coefficients of underlying 3-step explicit LM method are given by «о = 1,^1 = -1,^2 = «з = 0,^о = 0,^1 = 1, в2 = §, вз = —1- This method is of the second order and
zero-stable. However, its second characteristic polynomial is not stable. Finally, we carry out a similar experiment with well-known 2-step Adams-Moulton (AM2) method, which is of the third order and whose second characteristic polynomial is not stable. The numerical results are presented in Table 4. If it is necessary, for the approximation of retarded values, 4- and 5-node forward Newton interpolations can be used, respectively. It is clearly seen that HELM3 method and AM2 method for (6) are convergent of the second order and third order, respectively, but the direct discretizations based on the same underlying methods fail, which verifies the discussion in Remark 2.
Table 3
Numerical results for (21) on interval [0, 20] with A = -1,5,u = 10, a = 0, 5, b = 1,c = 0, 8
HELM3 method Standard HELM3 method
h = 0,1 Error in x1 Rate in xi Error in x1 Rate in x1
h 4,6970e-01 _ 6.6815c 1 15 ............
h/2 7,7009e-02 2,6086e 1.37!2c 270
h/4 l,6818e-02 2,1950 7.270 Ic 305 Not
h/8 4,1021e-03 2,0356 l,9288e+305 convergent
h/16 l,0138e-03 2,0165 2.9598c 30 !
h/32 2,5205e-04 2,0080 9.1235c 303
h = 0,1 Error in x2 Rate in x2 Error in x2 Rate in x2
h l,4985e-02 ............ 3.3256c 1 13
h/2 3,4649e-03 2,1126 6.8368c 267
h/4 8,3080e-04 2,0603 5,9109e+303 Not
h/8 2,0322e-04 2,0314 2.9392c 303 convergent
h/16 5,0236e-05 2,0163 8.5637c 302
h/32 l,2487e-05 2,0082 5.0161c 302
Table 4
Numerical results for (21) on interval [0, 20] with A = -1,5,u = 10, a = 0, 5, b = 1, c = 0, 8
AM2 method Standard AM2 method
h = 0,1 Error in x1 Rate in x1 Error in x1 Rate in x1
h 1.21 1 lc-03 ............ 1.71 11c 125 ............
h/2 l,4609e-04 3,0518 1.2302c 172
h/4 l,7941e-05 3,0256 3.5169c 266 Not
h/8 2.227 lc-06 3,0100 1.1311c 305 convergent
h/16 2,7735e-07 3,0054 Inf
h/32 3,4612e-08 3,0024 Inf
h = 0,1 Error in x2 Rate in x2 Error in x2 Rate in x2
h 5,9310e-05 ............ 8.5280c 122 ............
h/2 7,2103e-06 3,0402 6.1205c 169
h/4 8,8852e-07 3,0206 1.7616c 261 Not
h/8 l,1031e-07 3,0098 8.323 Ic 302 convergent
h/16 l,3741e-08 3,0050 Inf
h/32 l,7147e-09 3,0025 Inf
Conclusion
In this work, we have analyzed the conditioning of IVPs for a class of structured strangeness-free DDAEs with constant delay (3). Then, we have proposed LM methods combined with interpolation for solving this class of DDAEs. Convergence of the numerical methods have been established. The numerical experiments have confirmed the theoretical results.
There are numerous difficulties that arise in solving general DDAEs. From the analysis and results of this paper, it seems that there are many interesting works in the future. First, we should analyze the use of Runge-Kutta methods with continuous extension for DDAEs. Secondly, we could extend the analysis of numerical methods for general strangeness-free DDAEs with non-constant (time-varying, state-dependent) delay, which are much more difficult than the constant-delay case and would require more efforts.
Acknowledgements. M. V. Bulatov's work has been supported by the Russian Foundation for Basic Research, Grant Nos. 18-01-00643, 18-51-54001; V.H. Linh and N.D. Truong's work has been supported by Nafosted Project No. 101.02-2011.314-
References
1. Baker C.T.H., Paul C.A.H., Tian H. Differential Algebraic Equations with After-Effect. Journal of Computational and Applied Mathematics, 2002, vol. 140, pp. 63-80.
2. Bellen A., Maset S., Zennaro M., Guglielmi N. Recent Trends in the Numerical Solution of Retarded Functional Differential Equations. Acta Num.erica, 2009, vol. 18, pp. 1-110.
3. Ha P. Analysis and Numerical Solution of Delay Differential-Algebraic Equations. PhD Thesis. Berlin, 2015.
4. Shampine L.F., Gahinet P. Delay-Differential-Algebraic Equations in Control Theory. Applied Numerical Mathematics, 2006, vol. 56, pp. 574-588.
5. Bellen A., Zennaro M. Numerical Methods for Delay Differential Equations. Oxford, Oxford University Press, 2003.
6. Bellman R., Cooke K.L. Differential-Difference Equations. N.Y., Academic Press, 1963.
7. Ascher U., Petzold L. Computer Methods for Ordinary Differential Equations and Differential-Algebraic Equations. Philadelphia, SIAM Society for Industrial and Applied Mathematics, 1998.
8. Hairer E., Wanner G. Solving Ordinary Differential Equation II. Berlin, Springer-Verlag, 1996.
9. Kunkel P., Mehrmann V. Differential-Algebraic Equations Analysis and Numerical Solution. Zurich, EMS Publishing House, 2006.
10. Ha P., Mehrmann V., Steinbrecher A. Analysis of Linear Variable Coefficient Delay Differential-Algebraic Equations. Journal of Dynamics and Differential Equations, 2014, vol. 26, pp. 1-26.
11. Ascher U., Petzold L. The Numerical Solution of Delay-Differential-Algebraic Equations of Retarded and Neutral Type. SIAM Journal on Numerical Analysis, 1995, vol. 32, pp. 1635-1657.
12. Hauber R. Numerical Treatment of Retarded Differential-Algebraic Equations by Collocation Methods. Advances in Computational Mathematics, 1997, vol. 7, pp. 573-592.
13. Hongliang L., Aiguo X. Convergence of Linear Multistep Methods and One-Leg Methods for Index-2 Differential-Algebraic Equations with a Variable Delay. Advances in Applied Mathematics and Mechanics, 2012, vol. 4, no. 5, pp. 636-646.
14. Linh V.H., Mehrmann V. Efficient Integration of Matrix-Valued Non-Stiff DAEs by Half-Explicit Methods. Journal of Computational and Applied Mathematics, 2014, vol. 262, pp. 346-360.
15. Linh V.H., Truong N.D. Stable Numerical Solution for a Class of Structured Differential-Algebraic Equations by Linear Multistep Methods. Submitted for publication, 2018.
16. Bulatov M.V., Linh V.H., Solovarova L.S. On BDF-Based Multistep Schemes for Some Classes of Linear Differential-Algebraic Equations of Index at Most 2. Acta Mathematica Vietnamica, 2016, vol. 41, no. 4, pp. 715-730.
17. Linh V.H., Truong N.D. Runge-Kutta Methods Revisited for a Class of Structured Strangeness-Free DAEs. Electronic transactions on numerical analysis ETNA, 2018, vol. 48, pp.131-155.
18. Ascher U., Petzold L. Stability of Computational Methods for Constrained Dynamics Systems. SIAM Journal on Scientific Computing, 1993, vol. 14, pp. 95-120.
19. Hairer E., Norsett S.P., Wanner G. Solving Ordinary Differential Equations I. Berlin, Springer-Verlag, 1996.
Received August 08, 2018
УДК 519.62 Б01: 10.14529/ттр180406
АНАЛИЗ СХОДИМОСТИ ЛИНЕЙНЫХ МНОГОШАГОВЫХ МЕТОДОВ ДЛЯ РЕШЕНИЯ ОДНОГО КЛАССА ДИФФЕРЕНЦИАЛЬНО-АЛГЕБРАИЧЕСКИХ УРАВНЕНИЙ
Ву Хоанг Линь1, Нгуен Дуй Труонг2, М.В. Булатов3
Национальный университет Вьетнама, г. Ханой, Вьетнам 2Университет Тран Чок Туан, г. Ханой, Вьетнам
3Институт динамики систем и теории управления им. В.М. Матросова СО РАН, г. Иркутск, Российская Федерация
Дифференциально-алгебраические уравнения (ДАУ) с запаздываниями используются для моделирования реальных явлений, в которых могут одновременно присутствовать ограничения и запаздывания. Известно также, что решение ДАУ с запаздываниями является более сложной задачей, чем решение ДАУ без запаздываний, т.к. в случае с запаздываниями обычно требуется приближение решений на предыдущих временных отрезках и часто можно наблюдать разрыв у старших производных решений. В последнее время нами были предложены линейные многошаговые методы
решения для ДАУ низкого индекса без запаздывания. В данной работе мы расширили применение разработанных методов и используем их для решения ДАУ высокого индекса с постоянным запаздыванием. Для аппроксимации решений при запаздывании используется полиномиальная интерполяция. Представлен анализ сходимости линейных многошаговых методов. Показано, что, как и в случае отсутствия запаздывания, если вместо исходного ДАУ с запаздыванием мы дискретизируем особым образом переформулированное ДАУ, то для сходимости методов не требуется строгая устойчивость второго характеристического многочлена, поставленного в соответствие используемым методам. Теоретические выкладки проиллюстрированы численными расчетами.
Ключевые слова: дифференциально-алгебраические уравнения с запаздыванием; линейные многошаговые методы; устойчивость; сходимость.
Ву Хоанг Линь, профессор, Национальный университет Вьетнама (г. Ханой, Вьетнам), [email protected]. уп.
Нгуен Дуй Труонг, аспирант, Университет Тран Чок Туан (г. Ханой, Вьетнам), [email protected].
Михаил Валерьянович Булатов, доктор физико-математических наук, профессор, главный научный сотрудник, Институт динамики систем и теории управления им. В.М. Матросова СО РАН (г. Иркутск, Российская Федерация), [email protected].
Поступила в редакцию 8 августа 2018 г.
