- Email: [email protected]
Exact and discretized stability of the pantograph equation
APPLIED NUMERICAL MATHEMATICS ELSEVIER
Applied Numerical
Mathematics
24 ( 1997) 295-308
Exact and discretized stability of the pantograph equation * Arieh Iserles ’ Department of Applied Mathematics
and Theoretical Physics, University of Cambridge, Silver Street, Cambridge CB3 9EW, UK
1. Introduction The pantograph equation y’(t) = Ay(t)
+ BY(&) + CY’(@),
t 3 0,
(1.1)
Y(O) = Yo,
A, B, C are d x d complex matrices, ye E Cd and 4 E (0, l), possesses large number of applications in applied mathematics and engineering [14] and is fast becoming the standard model problem in the design of numerical methods for functional-differential equations with unbounded delays [Z-5,12,17,21]. There exists a natural inclination to extrapolate from the known theory of linear functionaldifferential equations with constant delay,
where
y’(t) = Ay(t) + By@ - T) + Cy’(t - T),
t 3 0,
(1.2)
where r > 0, and an initial value is given for all t E (-7, 01,to guess the behaviour of (1.1). This inclination should be resisted, because (1 .I) and (1.2) represent diverse mathematical structures and pose completely different mathematical challenges. On the one hand, the solution of (1.2) exhibits a familiar pattern of derivative discontinuities, while the solution of (l.l), when unique, is infinitely smooth. Moreover, the asymptotic stability condition for (1.2), namely that all the zeros of det [zI - A - ePTZ(B + ZC)] = 0,
.z E @,
reside in the left half-plane, is very different from the equivalent condition for the pantograph (cf. Section 2). Surprisingly enough, the constant-delay asymptotic stability condition is, at least in principle, more difficult to verify, although simplified conditions exist in special cases (cf., for example, [13]). Lest there is an impression that the constant-delay equation (1.2) is somehow more difficult than (1. l), we hasten to recount the other side of the story. While (1.2) is a relatively straightforward interval map, the pantograph (I. 1) is altogether more interesting: it is a d-dimensional dynamical system at * Based on a talk presented at the Volterra Centennial, the Second International Volterra and Delay Equations, Arizona State University, 27-30 May 1996. ’E-mail: [email protected]. 0168-9274/97/$17.00
0 1997 Elsevier Science B.V. All rights reserved.
PII SO168-9274(97)00027-5
Conference
on the Numerical
Solution of
296
A. Iserles /Applied
Numerical Mathematics
24 (1997) 295-308
t = 0 but an infinite-dimensional one for every t > 0. Its behaviour, being determined by a countable number of modes, is more akin to nonlinear differential equations. In particular, the equation displays a formidable richness of interesting phenomena on its stability boundary. The intrinsic complexity of the pantograph equation is further emphasized under discretization. While a standard finite-difference formula, whether multistep or Runge-Kutta, renders (1.2) into a stationary (albeit, high-degree) difference equation, the pantograph (1.1) yields a nonstationary system, whose analysis is significantly more difficult. Indeed, the results to date are incomplete and their derivation calls for new mathematical techniques. The theme of the paper is the range of phenomena which the pantograph equation-whether the exact or the discretized-displays inside and on its stability boundary. In Section 2 we describe the behaviour of the exact equation, while Section 3 is devoted to its discretization. We hasten to emphasize that, insofar as the numerical analysis of (1.1) is concerned, the results to date are incomplete and open problems abound. Indeed, one of the goals of this brief survey is to single out the many interesting open problems and challenges for the attention of the broader numerical community. The pantograph equation (1.1) is a fruitful point of departure for generalizations: (1) The case q > 1 has attracted significant attention because of its many applications in probability theory, theory of spin glasses and, perhaps most intriguingly, theory of wavelets [9,10]. (2) There is a growing measure of interest in q E @, because of its connection with self-similar dressing chains of Schrijdinger operators [ 11,21,22]. (3) Although only preliminary and tentative results are available with regard to nonlinear functionaldifferential equations with proportional delays, it is quite apparent that even the humble equation
y’(t) = QY@> +
Y(@)(l
-
t 3 0,
YW)>
Y(O)= Yo,
displays a whole range of intriguing phenomena [ 151. (4) It is natural to generalize (1.1) by considering equations with several proportional delays. As a matter of fact, it is perfectly viable to consider a continuum of delays-in other words, integro-differential equations of the form 1
1
y’(t)
= a~(4
+
Jy(G) d/&d + Jy’(qt)ddt),
t 3 0,
Y(O)
=
1,
0
0
where both dp and dv are complex-valued Bore1 measures. Such equations are interesting inter because they admit generalized hypergeometric solutions [ 181.
alia
2. The exact equation Let us suppose that a(C), the spectrum of the matrix C, does not contain points of the form q-’ for ! E Z+. Then the solution of (1.1) exists [ 141. Moreover, if p(C) < 1 then the solution is unique
WI. It is instructive equation y’(t)
=
to examine what happens if the above conditions
q/(t) +
by(@)
+
Y'(@),
t 2 0,
Y(O)
= 1,
fail. We first consider the scalar
(2.1)
A. lserles /Applied Numerical Mathematics 24 (1997) 295-308
291
which is typical of the general case que E o(C). It is easy to verify directly that a power series solution to (2.1) exists if and only if a + b = 0 and that in the latter case every function y(t) = (1 - r7) + 77eat, rl E C, is a solution. The equation y’(t) = a~/@>+ by(&) + cy’(&),
t 3 0,
Y(O) = 1,
(2.2)
possesses an infinite number of solutions whenever, a, b # 0 and /cl > showing that (2.2) admits two different Dirichlet expansions. We defer the until we have introduced this important concept, later in this section. We next proceed to analyze asymptotic stability. Considering the simple c = 0, we proceed by the classical technique of direct upper bounds, not since it exposes the shortcoming inherent in this approach. Thus, let
1. This can be proved by discussion of this equation scalar equation (2.2) with because it is effective but
Since i &ls@)/*
= Re[%t)y’@)]
= (Rea)]G)]*
+Re[@(MqQ]
G (Rea)l3+)]’
+ MN
it follows that iv’ < (Rea + Ibl)q. Therefore
q(t) < exp[2t(Rea
+ lbl)], t 3 0, and the solution tends to zero if
lb1< -Rea.
Rea < 0,
(2.3)
The main problem with (2.3) is that it constitutes but a small subset of the genuine domain of asymptotic stability. The asymptotic behaviour of the pantograph equation (1.1) has been analyzed, often in a simplified form, by a number of authors [7-9,14,19,21]. The most powerful tool in this endeavour is an expansion of the solution in Dirichlet series of the form y(t) = jFJ DneQ”tAv,
(2.4)
t > 0,
n=O
where each D, is a d x d complex matrix and v E Cd. To argue that (2.4) exists we stipulate that A is nonsingular and that X # qep for every A, I_LE a(A), and ! E N. Substitution of the expansion (2.4) into (1.1) results in linear recurrences AD, - qnD,A = --B&-l
- qn-lCD,_IA,
n E N.
Each step of the recursion entails a solution of a linear system and all such commutator equations are nonsingular because of the aforementioned condition on the eigenvalues of A. Setting Da = I and letting D(t) = C,“=. D,P, it has been demonstrated that det
D(t)
=
fi
det(I - tqnC)
n=Odet (I + tq”BA-l)
’
t a ”
298
A. Iserles /Applied
Numerical Mathematics
24 (1997) 295-308
therefore detD(0) = 1 and we let w = D-‘(0)ya. This proves the existence of the Dirichlet series (2.4) [14]. It is possible to prove that (2.4) converges for any t 3 0 if p(A-‘I?) < 1 (otherwise Dirichlet series need be replaced by Dirichlet-Taylor series [14]). Moreover, convergence is uniform for t E [0, cc) as long as all the eigenvalues of A reside in the complex left half-plane. However, the great virtue of the expansion (2.4) is that it tells the asymptotic behaviour of (1.1) at a glance. In particular, if Rea(A)
p(A-‘B)
< 0,
< 1
(2.5)
then limt,, y(t) = 0 and the zero solution is asymptotically stable [14]. Note that (2.5) is not always necessary for asymptotic stability. For example, the equation y’(t) = AY@) - CAY@>
+ CY'(@)
is solved by etAy(0), hence asymptotic stability requires just Rea(A) < 0. However, in a welldefined sense the condition (2.5) is “generically” necessary and sufficient for asymptotic stability (cf. the discussion in [21]). Specializing_(2.4) to Eq. (2.2), we derive for b # 0 the explicit expansion y(t)
=
neq”at,t > -go ( > ’ ’
(-UC/k Q)m O”(-4hdn
b
(w>n
a
hd2c where the Gauss-Heine
0
(2.6)
symbols (z; q)n are defined by the formula
n-l (ziq)n
=
n
(I
-
8z),
z,qEC,
nEz+U{cxl}.
k=O
In particular, the asymptotic Rea < 0,
stability condition
(2.5) becomes (2.7)
Ibl < Ial.
An interesting observation is that (2.7) is independent of c. An even more profound observation, though, is that the set (2.7) is usually much greater than the set (2.3). Fig. 1 displays all the values of b consistent with (2.5) (the inner disc) and with (2.7) (the outer disc) for a given value of a E @. Note that the nearer a is to the pure imaginary axis, the larger the discrepancy between the two sets. Before proceeding to examine the stability boundary, we discharge the earlier promise of considering the case ICI > 1. The function (2.6) remains a valid solution of (2.2), but it is possible to show that other solutions exist. Specifically, direct substitution and an elementary calculus with basic hypergeometric functions affirm that the Dirichlet expansion y(t)
=
(UC/kq)cc 2 (-Uclb;q)n (-w/b; c&x (4;dn n=O
0c 4
ne-q-nbt/~,t > H
0
>
(2.8)
is another solution of the same equation. Note that the expansion (2.8) converges for (cl > q but its derivative converges only for ICI > 1. Moreover, (2.8) is differentiable only a finite number of times-it can be proved that the only C” solution of (2.2) is the expansion (2.6) [21]. The boundary of the stability domain (2.5) consists of several parts. Firstly, in the case Rea(A) < 0, p(A-‘B) = 1, provided that all the eigenvalues of A-‘B with unit modulus have the same
299
A. Iserles /Applied Numerical Mathematics 24 (1997) 295-308
Imb
Reb Fig. I. Values of b E @ consistent dotted line denotes -Re a.
I
with (2.3) (the smaller disc) and (2.7) (the larger disc) for a given value of a E c. The
I
-5
5
0
x lob y’(t)
=
iy(t)
+
Fig. 2. Almost-periodic
geometric and algebraic exists t, > 0 such that I~?/@)- W)
solutions of the pantograph
multiplicity,
(/ < E,
y’(t)
f$jiy(t/lO)
=
iy(t)
+
equation (2.2) in the
there exists a periodic
function
$$jy(9t/lO)
(Rey, Imy) plane.
C$ and for every E > 0 there
t 3 tE.
Otherwise the solution of (1.1) is unbounded
in [0, co) for some initial values.
300
A. Iserles /Applied Numerical Mathematics 24 (1997) 295-308
-1
0
0
I
y’(t) = iy(t) - $y’(t/7)
7
y’(t) = iy(t) - aY’(2tlll)
I
1
0
0
.I
-1
y’(t) = iy(t) - $y(t/4)
- $iy’(t/4)
y’(t) = iy(t) - +y(2t/5)
- $iy’(2t/5)
Fig. 3. Rotational almost-symmetry in the exact solution of the pantograph equation (2.2). Note the dependence
on the value
of k!- k, where q = k/e.
The situation is, if at all, even more interesting max { Re A: X E a(A)}
when
= 0.
If det A = 0 then the solution is unbounded in [0, oc) for some initial values. Otherwise there exists an almost-periodic function 1c, and for every E > 0 there exists s, > 0 such that [(z!(t) - N>ll
< E,
t > SE.
We remind the reader that a function ~6 > 0 such that (($J@ + Q) - $Q)((
< 4
+ E L,[O, CG) is almost-periodic
if for every S there exists
t 3 0,
and liminfalc v6 > 0. If all the eigenvalues of A reside on ilk! then the solution y itself is almost-periodic. for example, for Eq. (2.2) when a E ill%\(O). Th e solutions of two such equations
This is the case, are displayed in
A. Iserles /Applied Numerical Mathematics 24 (1997) 295-308
301
Fig. 2 in the plane (Re y(/,Im y). Both are suggestive of a great deal of structure. Indeed, the little we know to date on the dynamics of almost-periodic functions of this form is certainly quite interesting. If 4 = k/e, where Ic, e E N are relatively prime, 1 < k 6 t - 1, then y possesses an almost-symmetry (a concept that parallels almost-periodicity) under rotations by integer multiples of 2~/([ - k) [ 141. This can be seen quite vividly in Fig. 3, which displays four solutions for different rational values of q. Further features of the solution can be discerned in numerical experimentation but, as things stand, they lack complete verification. Thus, again for q = k/t, where k and f? are relatively prime, each plot is composed of ! “strands”. Each such “strand” becomes itself, upon zooming, a braid of ! “strands” and so on. If q irrational the solution seems to be space-filling in an annulus [14].
3. The numerical solution It is relatively easy to extend standard ODE methods to the pantograph equation (1.1). The only extra ingredient is the provision of continuous output but this can be always ensured by interpolation or by using natural Runge-Kutta methods. In a sense, it is easier to design time-stepping algorithms for (1.1) than for (1.2), because the latter equation exhibits derivative discontinuities. Having said this, stability analysis of the discretized pantograph presents a formidable challenge. We mention in passing that perhaps the most efficient means of discretizing the pantograph equation is provided by a Dirichlet expansion. Yet, as long as we consider (1.1) as a model equation and wish to infer from its analysis to the broader issue of the numerical analysis of functional differential equations, we are not allowed to exploit features which are characteristic of the pantograph. The familiar trapezoidal rule is a good starting point for the analysis of standard numerical methods, although Runge-Kutta schemes present an added complexity [6]. Interpolating g(qt) by a piecewiselinear function on an equispaced grid results in the nonstationary recurrence relation
(3.1) where c2,1=q(n+i)-LqnJ,
nEZ+,
which approximates (2.2) with a local error of O(h3). Asymptotic behaviour of the method (3.1) can be analyzed by classical methods, but this approach suffers all the shortcomings of similarly naive estimates for the exact solution of (2.2) that have been already debated in Section 2 and that have led to the set (2.3) [5]. In a special case when q is a reciprocal of an integer, though, different analysis has been presented by Martin Buhmann and the author and it leads to a result consistent with both numerical experiment and intuition [4]. Cf. [20] for an extension of this analysis to Q-methods.
302
A. Iserles /Applied Numerical Mathematics 24 (1997) 295-308
In the case q = $ the recurrence Y2n+l = By2n
(3.1) reads
$5(3%X + Yn+1) + Z(Yn+l
+
Y2n+2 = SY2nfl
+ 3y,+1)
+ $(Yn
+ Z(yn+l
- Yn),
n E Z+,
(3.2)
- yn),
where B=
1+ kha
6=
1 - iha’
hb
x=
1 - iha’
c 1 - kha’
Although (3.2) can be analyzed in the above form, we turn our attention to Y2n+l = ~Y~TI + &yn,
n E Z+,
(3.3)
y2?2+2 = XYzTz+1 + P(Yn + Yn+1),
where X, p E @. Eq. (3.3) is a convenient model for difference equations with proportional delay. The investigation of its stability boundary is somewhat easier (and certainly simpler to present) than the corresponding theory of (3.2). Yet, to the extent of our present (and incomplete) understanding, the two equations share similar qualitative features on their stability boundaries. Intuitively speaking, X should be identified with % and ,u with ;S. Letting xn := yn+i - ynln,we obtain from (3.3) the recurrence Zo=X+2~-1, x2,+1
= XZ2n + wn
x2,+2
= XX2nSl
n E Z+.
+ w,
1
Let us denote by Z?(e), 0 E [0,27r], th e F ourier trunsfomz of the sequence {x,},,=~+. to prove that Z obeys the functional equation
E(8) = l _‘x,
[x:0 + ~(1 + eie)Z(20)] (mod27r),
0 E [0,2-/r].
It is possible
(3.4)
The function Z?is exceedingly complicated even for “nice” values of X and p, as should be evident from Fig. 4. The situation is even worse for (X( = 1, since then the function is unbounded. Yet, it is not outside the realm of analysis! Applying (3.4) repeatedly, we expand the Fourier transform into infinite series, q(j)
=
(A
+
Q
_
1)
5 l-IT=0’ u+e2kie) m, f9 E P,27d, m=O
nTzo( 1 - Xe2kie) ’
and use the mean ergodic theorem to deduce that Z? is an L2 [0,27r] function for all 1x1 < 1 and 21~1 < 11 - X(. From this there is but a short step to argue that also the Fourier transform of {Y~}~~z+ is L2[0,27r], whence yin + 0 and the recurrence (3.3) is asymptotically stable. Similar analysis, albeit more technically convoluted, applies to the numerical method (3.2), whereby the asymptotic stability condition is I??% < 1, 16:) < 11 - 81, which translates exactly into the condition (2.5). Hence, at least when q is a reciprocal of an integer, the trapezoidal rule recovers the correct asymptotic stability domain of the pantograph equation [4]. (Note that in the neutral case 0 < (cl < 1 we need to restrict the step-size so that ]c f ihbl < 11 - ihal.)
303
A. Iserles /Applied Numerical Mathematics 24 (1997) 295-308
-5
1
2
3
4
5
6
1
2
3
4
5
6
0
20
I
10
0
1'
- 10
0
1
2
3
4
5
6
0
1
2
3
4
5
6
Fig. 4. Real and imaginary parts of the Fourier transform ?? corresponding
to different choices of X and p.
The behaviour of (3.3) on the stability boundary 1x1 = 1 is considerably more complicated. However, a number of results are available and they indicate that the attractor of the sequence {Y~,}~~~+ can be an exceedingly complicated object, as can be seen in Fig. 5. Unfortunately, virtually all (bar numerical experiments) that is known in precise mathematical terms with regard to the attractor refers to the special case X = - 1, which does not correspond to a realistic numerical method (since 9%= 1 + ha + 0(h2)). We expand each yn in powers of p, yn = C-V
+ ~a,,&, e=i
n E Z+.
304
A. lserles /Applied
Numerical Mathematics
24 (1997) 29.5-308
r
c
3 I L
X = -1, p
zz fi
X=i,p=$
Fig. 5. Attractors of the soldon sequence of (3.3) for different values of X and ,u, sketched in the plane (Re yn, Im yn) Let
US commence
y1 = -1
+2/i,
y2 = +1
-2p
by examining
the first few values
+2/_L2,
y3 =-I
++*,
y4 =+I
-2/L*
+2/x3,
-2/L*
+2p3,
y5 = -1
+2/i
35 = +I
-2/L
+2p3,
y7=-1
+G3
>
y/8 =+I
-2p3
+ 2y4,
y9 = -1
+2/l
-2p3 + 2p4,
y10 = +1
-2p
+2/L* -2p3 + 2p4 >
Y/II = -I
+2p2
-2p3 + 2p4,
Y12 = +I
-2/L*
-I- 2p?.
of pn,
A. Iserles /Applied Numerical Mathematics 24 (1997) 29.5-308
305
Obviously, we can discern a pattern! Firstly, it appears that an,e E (0, +2, -2) and, secondly, there is clearly a binary structure to the sequence. It is not difficult to show that this is indeed the case and that the explicit form of the coefficients is
0, “2e+Im+j.e
=
j =0,1,...,2e-l j=2e-l,...,2’j =2$...>3~2~-‘j=3x2”-‘....>2e+‘-l
2, -2, 0:
- 1: 1. 1.
for all e E N and m E Z+ [16]. Let A, be the w-set of the sequence {Yn}nEZ+. In other words, A, is the closure of all the accumulation points of the sequence. The explicit form of the coefficients allows to evaluate specific members of A,. For example, as long as jl_~J< 1 it is true that lim Y(CI)/X T.403
= -G
1-P
and
We present next an explicit evaluation 2’-1=2@l,+j
J&y2+1)/3
=
1-P
1 +pL’
of yzT_ 1. Since for all e = 1,2, . . . , T - 1 it is true that
wherem=2r-e-1-l,
j12~~l-1,
we deduce that ~~_l,e = 0 in that range. Moreover, for ! = T we have m = 0, j = 2’ - 1, whence CQ _ i,r = 2, and trivially Q_ 1,e = 0 for 1 3 T + 1. Therefore ~2~1 = - 1 + 2~’ and lim ~7-1
r+x
= -1.
It is possible to prove that, in general, each A, can be inscribed into a complex disc of radius (1 - r_ll/(l - 1~~1)[16]. It is known that A0 = {-l,+l} and A l/2 = [- 1, 11. Perhaps more surprising result is that every element of the sequence {Y%}~,=~+is in Ap. In order to prove this we pick an arbitrary n E Z+ and let rnk := n+2 Its-r, k E Z+. Obviously, the mk’s form a pairwise-distinct sequence and, provided that k is large enough, binary expansions of n and rnk differ only in the kth digit. It is easy to prove that arnlc,e = an,e, ! = 0, 1. . . ) k - 1, whence IYn-YlnkI s25p,~=3po B=k
for all Ir_~l< 1 [16]. As a matter of fact, unless CL= 0, the set A, is of an infinite cardinality. Moreover, it follows at once from the definition that there exists a countable number of components of A, in every nontrivial neighbourhood of every element of this set. Fig. 6 displays the sets A, for four choices of p E C. It is indicative of fractal sets, as is the above analysis. We proceed next to identify A,, p # 0, with a different mathematical construct, thereby proving that our suspicion is well-founded: the attractors indeed possess a fractional Hausdorff dimension. Letting X = -1 in (3.3), we eliminate pyn and pyn+i from these expressions. The first equation can be rewritten in the form ;PL(Y2n
PLY72 =
and substitution Y2nfl
-
+ Y2n+l)
in the second equation yields
Y2n =
y2n-1
-
Y2n-2.
306
A. Iserles /Applied
Numerical Mathematics
24 (1997) 295-308
I
0
0 p = +(l$i)
$i
,u =
1
1
0 -1
-10
0i? 4 $
I
J
0
p=
p=$-ai
&(l+
i)
Fig. 6. w-sets A, for different values of p E &1, displayed in the plane (Red,,
Consequently,
Imd,).
by induction,
Y2n+l - Y2n = Yl - Yo = 41
- p).
Adding and subtracting the first equation results in Y2n = PYn + (1 - P),
n E N,
Y2nSl
n E z+,
= PYn -
respectively. fl(4
(1 - P),
(3.5)
Letting := W + (1 - 4,
the recurrence
f2(4
:=
llz
-
(1
-
4,
z E c,
(3.5) becomes
Y2n = fl (Yn>,
Y2nfl
= f2(Yn).
This, in turn, has the same attractor as a dynamical system which is obtained when for each n E Z4 we obtain yn+t from yn by iterating either ft or f2 with equal probability [l]. The attractor .A, is thus identified with a probabilistic mixture of Julia sets, whose properties have been a subject to a thoroughgoing mathematical analysis in the last decade.
A. Iserles /Applied Numerical Mathematics 24 (1997) 295-308
307
To state that the present body of knowledge on the discretized pantograph equation is incomplete constitutes an understatement. Insofar as asymptotic stability is concerned, the method of proof extends just to Q which is a reciprocal of an integer and on the stability boundary we can explain only the nonrealistic case X = - 1. Yet, computer experimentation indicates that our results represent general patterns. It is fair to formulate the conjecture that, subject to further step-size restrictions in the neutral case, the discretization of (2.2) by an arbitrary multistep method is stable whenever (bl < la( and ha belongs to the linear stability domain of the underlying ODE scheme. Moreover, computer experiments seem to indicate that values of ha on the boundary of the linear stability domain lead to fractal attractors. Only limitations of our current methods of proof and of our understanding stand between us and a verification (or refutation) of these conjectures.
References [l] M.F. Barnsley and S.G. Den&o, Rational approximation of fractals, in: P.R. Graves-Morris, E.B. Saff and R.S. Varga, eds., Rational Approximation and Interpolation, Lecture Notes in Mathematics 1105 (Springer, Berlin, 1984) 73-88. [2] A. Bellen, N. Guglielmi and L. Torelli, Asymptotic stability properties of theta-methods for the pantograph equation, Quaderno Matematico 370 de1 Dipartimento di Scienze Matematiche dell’Universit8 di Trieste (1996). [3] H. Brunner, On the discretization of differential and Volterra integral equations with variable delay, Technical Report, Memorial University of Newfoundland (1996). [4] M.D. Buhmann and A. Iserles, On the dynamics of a discretized neutral equation, ZMA J. Numer. Anal. 12 (1992) 339-363. [5] M.D. Buhmann and A. Iserles, Stability of the discretized pantograph differential equation, Math. Comp. 60 (1993) 575-589. [6] M.D. Buhmann, A. Iserles and S.P. Norsett, Runge-Kutta methods for neutral differential equations, in: R.P. Agarwal, ed., Contributions in Numerical Mathematics, World Scientific Series in Applicable Analysis (World Scientific, Singapore, 1993) 85-98. [7] J. Carr and J. Dyson, The functional differential equation y’(x) = ay(Xz)+by(z), Proc. Roy. Sot. Edinburgh Sect. A 74 (1974/1975) 165-174. [8] J. Carr and J. Dyson, The matrix differential equation y’(x) = Ay(Xz) +By(z), Proc. Roy. Sot. Edinburgh Sect. A 75 (1975/1976) 5-22. of solutions of functional and differential-functional equations with several [91 G. Derfel, Behaviour transformations of the independent variable, Ukrai: Mat. Zh. 34 (1982) 350-356. [lo] G. Derfel, N. Dyn and D. Levin, Generalized refinement equation and subdivision processes, J. Approx. Theory 80 (1995) 100-111. [l l] G. Derfel and A. Iserles, The pantograph equation in the complex plane, DAMTP Technical Report 1996/NAOl (1996); also: J. Math. Anal. Appl., submitted. [12] L. Fox, D.F. Mayers, J.R. Ockendon and A.B. Tayler, On a functional-differential equation, J. Inst. Math. Appl. 8 (1971) 271-307. [13] K.J. in ‘t Hout, The stability of f&methods for systems of delay differential equations, Ann. Numer. Math. 1 (1994) 323-334. [ 141 A. Iserles, On the generalized pantograph functionaldifferential equation, European .I. Appl. Math. 4 (1993) l-38. [15] A. Iserles, On nonlinear delay differential equations, Trans. Amel: Math. Sot. 344 (1994) 441-477.
308
A. Iserles /Applied Numerical Mathematics 24 (1997) 295-308
[ 161 A. Iserles, The asymptotic behaviour of certain difference equations with proportional delay, Comput. Math. Appl. 28 (1994) 141-152. [17] A. Iserles and Y. Liu, On neutral functional-differential equations with proportional delays, DAMTP Technical Report 1993/NA3 (1993); also: J. Math. Anal. Appl., to appear. [18] A. Iserles and Y. Liu, On pantograph integro-differential equations, J. Integral Equations Appl. 6 (1994) 213-237. [19] T. Kato and J.B. McLeod, The functional-differential equation y’(z) = ay(Xx) + by(z), Bull. Amex Math. Sot. 77 (1971) 891-937. [20] Y. Liu, Stability of Q-methods for neutral functional-differential equations, Numel: Math. 70 (1995) 473485. [21] Y. Liu, On functional differential equations with proportional delays, Ph.D. Dissertation, University of
Cambridge ( 1996). [22] S. Skorik and V. Spiridonov, Self-similar potentials and the q-oscillator algebra at roots of unity, iAt. Math. Phys. 28 (1993) 59-74.
Applied Numerical
Mathematics
24 ( 1997) 295-308
Exact and discretized stability of the pantograph equation * Arieh Iserles ’ Department of Applied Mathematics
and Theoretical Physics, University of Cambridge, Silver Street, Cambridge CB3 9EW, UK
1. Introduction The pantograph equation y’(t) = Ay(t)
+ BY(&) + CY’(@),
t 3 0,
(1.1)
Y(O) = Yo,
A, B, C are d x d complex matrices, ye E Cd and 4 E (0, l), possesses large number of applications in applied mathematics and engineering [14] and is fast becoming the standard model problem in the design of numerical methods for functional-differential equations with unbounded delays [Z-5,12,17,21]. There exists a natural inclination to extrapolate from the known theory of linear functionaldifferential equations with constant delay,
where
y’(t) = Ay(t) + By@ - T) + Cy’(t - T),
t 3 0,
(1.2)
where r > 0, and an initial value is given for all t E (-7, 01,to guess the behaviour of (1.1). This inclination should be resisted, because (1 .I) and (1.2) represent diverse mathematical structures and pose completely different mathematical challenges. On the one hand, the solution of (1.2) exhibits a familiar pattern of derivative discontinuities, while the solution of (l.l), when unique, is infinitely smooth. Moreover, the asymptotic stability condition for (1.2), namely that all the zeros of det [zI - A - ePTZ(B + ZC)] = 0,
.z E @,
reside in the left half-plane, is very different from the equivalent condition for the pantograph (cf. Section 2). Surprisingly enough, the constant-delay asymptotic stability condition is, at least in principle, more difficult to verify, although simplified conditions exist in special cases (cf., for example, [13]). Lest there is an impression that the constant-delay equation (1.2) is somehow more difficult than (1. l), we hasten to recount the other side of the story. While (1.2) is a relatively straightforward interval map, the pantograph (I. 1) is altogether more interesting: it is a d-dimensional dynamical system at * Based on a talk presented at the Volterra Centennial, the Second International Volterra and Delay Equations, Arizona State University, 27-30 May 1996. ’E-mail: [email protected]. 0168-9274/97/$17.00
0 1997 Elsevier Science B.V. All rights reserved.
PII SO168-9274(97)00027-5
Conference
on the Numerical
Solution of
296
A. Iserles /Applied
Numerical Mathematics
24 (1997) 295-308
t = 0 but an infinite-dimensional one for every t > 0. Its behaviour, being determined by a countable number of modes, is more akin to nonlinear differential equations. In particular, the equation displays a formidable richness of interesting phenomena on its stability boundary. The intrinsic complexity of the pantograph equation is further emphasized under discretization. While a standard finite-difference formula, whether multistep or Runge-Kutta, renders (1.2) into a stationary (albeit, high-degree) difference equation, the pantograph (1.1) yields a nonstationary system, whose analysis is significantly more difficult. Indeed, the results to date are incomplete and their derivation calls for new mathematical techniques. The theme of the paper is the range of phenomena which the pantograph equation-whether the exact or the discretized-displays inside and on its stability boundary. In Section 2 we describe the behaviour of the exact equation, while Section 3 is devoted to its discretization. We hasten to emphasize that, insofar as the numerical analysis of (1.1) is concerned, the results to date are incomplete and open problems abound. Indeed, one of the goals of this brief survey is to single out the many interesting open problems and challenges for the attention of the broader numerical community. The pantograph equation (1.1) is a fruitful point of departure for generalizations: (1) The case q > 1 has attracted significant attention because of its many applications in probability theory, theory of spin glasses and, perhaps most intriguingly, theory of wavelets [9,10]. (2) There is a growing measure of interest in q E @, because of its connection with self-similar dressing chains of Schrijdinger operators [ 11,21,22]. (3) Although only preliminary and tentative results are available with regard to nonlinear functionaldifferential equations with proportional delays, it is quite apparent that even the humble equation
y’(t) = QY@> +
Y(@)(l
-
t 3 0,
YW)>
Y(O)= Yo,
displays a whole range of intriguing phenomena [ 151. (4) It is natural to generalize (1.1) by considering equations with several proportional delays. As a matter of fact, it is perfectly viable to consider a continuum of delays-in other words, integro-differential equations of the form 1
1
y’(t)
= a~(4
+
Jy(G) d/&d + Jy’(qt)ddt),
t 3 0,
Y(O)
=
1,
0
0
where both dp and dv are complex-valued Bore1 measures. Such equations are interesting inter because they admit generalized hypergeometric solutions [ 181.
alia
2. The exact equation Let us suppose that a(C), the spectrum of the matrix C, does not contain points of the form q-’ for ! E Z+. Then the solution of (1.1) exists [ 141. Moreover, if p(C) < 1 then the solution is unique
WI. It is instructive equation y’(t)
=
to examine what happens if the above conditions
q/(t) +
by(@)
+
Y'(@),
t 2 0,
Y(O)
= 1,
fail. We first consider the scalar
(2.1)
A. lserles /Applied Numerical Mathematics 24 (1997) 295-308
291
which is typical of the general case que E o(C). It is easy to verify directly that a power series solution to (2.1) exists if and only if a + b = 0 and that in the latter case every function y(t) = (1 - r7) + 77eat, rl E C, is a solution. The equation y’(t) = a~/@>+ by(&) + cy’(&),
t 3 0,
Y(O) = 1,
(2.2)
possesses an infinite number of solutions whenever, a, b # 0 and /cl > showing that (2.2) admits two different Dirichlet expansions. We defer the until we have introduced this important concept, later in this section. We next proceed to analyze asymptotic stability. Considering the simple c = 0, we proceed by the classical technique of direct upper bounds, not since it exposes the shortcoming inherent in this approach. Thus, let
1. This can be proved by discussion of this equation scalar equation (2.2) with because it is effective but
Since i &ls@)/*
= Re[%t)y’@)]
= (Rea)]G)]*
+Re[@(MqQ]
G (Rea)l3+)]’
+ MN
it follows that iv’ < (Rea + Ibl)q. Therefore
q(t) < exp[2t(Rea
+ lbl)], t 3 0, and the solution tends to zero if
lb1< -Rea.
Rea < 0,
(2.3)
The main problem with (2.3) is that it constitutes but a small subset of the genuine domain of asymptotic stability. The asymptotic behaviour of the pantograph equation (1.1) has been analyzed, often in a simplified form, by a number of authors [7-9,14,19,21]. The most powerful tool in this endeavour is an expansion of the solution in Dirichlet series of the form y(t) = jFJ DneQ”tAv,
(2.4)
t > 0,
n=O
where each D, is a d x d complex matrix and v E Cd. To argue that (2.4) exists we stipulate that A is nonsingular and that X # qep for every A, I_LE a(A), and ! E N. Substitution of the expansion (2.4) into (1.1) results in linear recurrences AD, - qnD,A = --B&-l
- qn-lCD,_IA,
n E N.
Each step of the recursion entails a solution of a linear system and all such commutator equations are nonsingular because of the aforementioned condition on the eigenvalues of A. Setting Da = I and letting D(t) = C,“=. D,P, it has been demonstrated that det
D(t)
=
fi
det(I - tqnC)
n=Odet (I + tq”BA-l)
’
t a ”
298
A. Iserles /Applied
Numerical Mathematics
24 (1997) 295-308
therefore detD(0) = 1 and we let w = D-‘(0)ya. This proves the existence of the Dirichlet series (2.4) [14]. It is possible to prove that (2.4) converges for any t 3 0 if p(A-‘I?) < 1 (otherwise Dirichlet series need be replaced by Dirichlet-Taylor series [14]). Moreover, convergence is uniform for t E [0, cc) as long as all the eigenvalues of A reside in the complex left half-plane. However, the great virtue of the expansion (2.4) is that it tells the asymptotic behaviour of (1.1) at a glance. In particular, if Rea(A)
p(A-‘B)
< 0,
< 1
(2.5)
then limt,, y(t) = 0 and the zero solution is asymptotically stable [14]. Note that (2.5) is not always necessary for asymptotic stability. For example, the equation y’(t) = AY@) - CAY@>
+ CY'(@)
is solved by etAy(0), hence asymptotic stability requires just Rea(A) < 0. However, in a welldefined sense the condition (2.5) is “generically” necessary and sufficient for asymptotic stability (cf. the discussion in [21]). Specializing_(2.4) to Eq. (2.2), we derive for b # 0 the explicit expansion y(t)
=
neq”at,t > -go ( > ’ ’
(-UC/k Q)m O”(-4hdn
b
(w>n
a
hd2c where the Gauss-Heine
0
(2.6)
symbols (z; q)n are defined by the formula
n-l (ziq)n
=
n
(I
-
8z),
z,qEC,
nEz+U{cxl}.
k=O
In particular, the asymptotic Rea < 0,
stability condition
(2.5) becomes (2.7)
Ibl < Ial.
An interesting observation is that (2.7) is independent of c. An even more profound observation, though, is that the set (2.7) is usually much greater than the set (2.3). Fig. 1 displays all the values of b consistent with (2.5) (the inner disc) and with (2.7) (the outer disc) for a given value of a E @. Note that the nearer a is to the pure imaginary axis, the larger the discrepancy between the two sets. Before proceeding to examine the stability boundary, we discharge the earlier promise of considering the case ICI > 1. The function (2.6) remains a valid solution of (2.2), but it is possible to show that other solutions exist. Specifically, direct substitution and an elementary calculus with basic hypergeometric functions affirm that the Dirichlet expansion y(t)
=
(UC/kq)cc 2 (-Uclb;q)n (-w/b; c&x (4;dn n=O
0c 4
ne-q-nbt/~,t > H
0
>
(2.8)
is another solution of the same equation. Note that the expansion (2.8) converges for (cl > q but its derivative converges only for ICI > 1. Moreover, (2.8) is differentiable only a finite number of times-it can be proved that the only C” solution of (2.2) is the expansion (2.6) [21]. The boundary of the stability domain (2.5) consists of several parts. Firstly, in the case Rea(A) < 0, p(A-‘B) = 1, provided that all the eigenvalues of A-‘B with unit modulus have the same
299
A. Iserles /Applied Numerical Mathematics 24 (1997) 295-308
Imb
Reb Fig. I. Values of b E @ consistent dotted line denotes -Re a.
I
with (2.3) (the smaller disc) and (2.7) (the larger disc) for a given value of a E c. The
I
-5
5
0
x lob y’(t)
=
iy(t)
+
Fig. 2. Almost-periodic
geometric and algebraic exists t, > 0 such that I~?/@)- W)
solutions of the pantograph
multiplicity,
(/ < E,
y’(t)
f$jiy(t/lO)
=
iy(t)
+
equation (2.2) in the
there exists a periodic
function
$$jy(9t/lO)
(Rey, Imy) plane.
C$ and for every E > 0 there
t 3 tE.
Otherwise the solution of (1.1) is unbounded
in [0, co) for some initial values.
300
A. Iserles /Applied Numerical Mathematics 24 (1997) 295-308
-1
0
0
I
y’(t) = iy(t) - $y’(t/7)
7
y’(t) = iy(t) - aY’(2tlll)
I
1
0
0
.I
-1
y’(t) = iy(t) - $y(t/4)
- $iy’(t/4)
y’(t) = iy(t) - +y(2t/5)
- $iy’(2t/5)
Fig. 3. Rotational almost-symmetry in the exact solution of the pantograph equation (2.2). Note the dependence
on the value
of k!- k, where q = k/e.
The situation is, if at all, even more interesting max { Re A: X E a(A)}
when
= 0.
If det A = 0 then the solution is unbounded in [0, oc) for some initial values. Otherwise there exists an almost-periodic function 1c, and for every E > 0 there exists s, > 0 such that [(z!(t) - N>ll
< E,
t > SE.
We remind the reader that a function ~6 > 0 such that (($J@ + Q) - $Q)((
< 4
+ E L,[O, CG) is almost-periodic
if for every S there exists
t 3 0,
and liminfalc v6 > 0. If all the eigenvalues of A reside on ilk! then the solution y itself is almost-periodic. for example, for Eq. (2.2) when a E ill%\(O). Th e solutions of two such equations
This is the case, are displayed in
A. Iserles /Applied Numerical Mathematics 24 (1997) 295-308
301
Fig. 2 in the plane (Re y(/,Im y). Both are suggestive of a great deal of structure. Indeed, the little we know to date on the dynamics of almost-periodic functions of this form is certainly quite interesting. If 4 = k/e, where Ic, e E N are relatively prime, 1 < k 6 t - 1, then y possesses an almost-symmetry (a concept that parallels almost-periodicity) under rotations by integer multiples of 2~/([ - k) [ 141. This can be seen quite vividly in Fig. 3, which displays four solutions for different rational values of q. Further features of the solution can be discerned in numerical experimentation but, as things stand, they lack complete verification. Thus, again for q = k/t, where k and f? are relatively prime, each plot is composed of ! “strands”. Each such “strand” becomes itself, upon zooming, a braid of ! “strands” and so on. If q irrational the solution seems to be space-filling in an annulus [14].
3. The numerical solution It is relatively easy to extend standard ODE methods to the pantograph equation (1.1). The only extra ingredient is the provision of continuous output but this can be always ensured by interpolation or by using natural Runge-Kutta methods. In a sense, it is easier to design time-stepping algorithms for (1.1) than for (1.2), because the latter equation exhibits derivative discontinuities. Having said this, stability analysis of the discretized pantograph presents a formidable challenge. We mention in passing that perhaps the most efficient means of discretizing the pantograph equation is provided by a Dirichlet expansion. Yet, as long as we consider (1.1) as a model equation and wish to infer from its analysis to the broader issue of the numerical analysis of functional differential equations, we are not allowed to exploit features which are characteristic of the pantograph. The familiar trapezoidal rule is a good starting point for the analysis of standard numerical methods, although Runge-Kutta schemes present an added complexity [6]. Interpolating g(qt) by a piecewiselinear function on an equispaced grid results in the nonstationary recurrence relation
(3.1) where c2,1=q(n+i)-LqnJ,
nEZ+,
which approximates (2.2) with a local error of O(h3). Asymptotic behaviour of the method (3.1) can be analyzed by classical methods, but this approach suffers all the shortcomings of similarly naive estimates for the exact solution of (2.2) that have been already debated in Section 2 and that have led to the set (2.3) [5]. In a special case when q is a reciprocal of an integer, though, different analysis has been presented by Martin Buhmann and the author and it leads to a result consistent with both numerical experiment and intuition [4]. Cf. [20] for an extension of this analysis to Q-methods.
302
A. Iserles /Applied Numerical Mathematics 24 (1997) 295-308
In the case q = $ the recurrence Y2n+l = By2n
(3.1) reads
$5(3%X + Yn+1) + Z(Yn+l
+
Y2n+2 = SY2nfl
+ 3y,+1)
+ $(Yn
+ Z(yn+l
- Yn),
n E Z+,
(3.2)
- yn),
where B=
1+ kha
6=
1 - iha’
hb
x=
1 - iha’
c 1 - kha’
Although (3.2) can be analyzed in the above form, we turn our attention to Y2n+l = ~Y~TI + &yn,
n E Z+,
(3.3)
y2?2+2 = XYzTz+1 + P(Yn + Yn+1),
where X, p E @. Eq. (3.3) is a convenient model for difference equations with proportional delay. The investigation of its stability boundary is somewhat easier (and certainly simpler to present) than the corresponding theory of (3.2). Yet, to the extent of our present (and incomplete) understanding, the two equations share similar qualitative features on their stability boundaries. Intuitively speaking, X should be identified with % and ,u with ;S. Letting xn := yn+i - ynln,we obtain from (3.3) the recurrence Zo=X+2~-1, x2,+1
= XZ2n + wn
x2,+2
= XX2nSl
n E Z+.
+ w,
1
Let us denote by Z?(e), 0 E [0,27r], th e F ourier trunsfomz of the sequence {x,},,=~+. to prove that Z obeys the functional equation
E(8) = l _‘x,
[x:0 + ~(1 + eie)Z(20)] (mod27r),
0 E [0,2-/r].
It is possible
(3.4)
The function Z?is exceedingly complicated even for “nice” values of X and p, as should be evident from Fig. 4. The situation is even worse for (X( = 1, since then the function is unbounded. Yet, it is not outside the realm of analysis! Applying (3.4) repeatedly, we expand the Fourier transform into infinite series, q(j)
=
(A
+
Q
_
1)
5 l-IT=0’ u+e2kie) m, f9 E P,27d, m=O
nTzo( 1 - Xe2kie) ’
and use the mean ergodic theorem to deduce that Z? is an L2 [0,27r] function for all 1x1 < 1 and 21~1 < 11 - X(. From this there is but a short step to argue that also the Fourier transform of {Y~}~~z+ is L2[0,27r], whence yin + 0 and the recurrence (3.3) is asymptotically stable. Similar analysis, albeit more technically convoluted, applies to the numerical method (3.2), whereby the asymptotic stability condition is I??% < 1, 16:) < 11 - 81, which translates exactly into the condition (2.5). Hence, at least when q is a reciprocal of an integer, the trapezoidal rule recovers the correct asymptotic stability domain of the pantograph equation [4]. (Note that in the neutral case 0 < (cl < 1 we need to restrict the step-size so that ]c f ihbl < 11 - ihal.)
303
A. Iserles /Applied Numerical Mathematics 24 (1997) 295-308
-5
1
2
3
4
5
6
1
2
3
4
5
6
0
20
I
10
0
1'
- 10
0
1
2
3
4
5
6
0
1
2
3
4
5
6
Fig. 4. Real and imaginary parts of the Fourier transform ?? corresponding
to different choices of X and p.
The behaviour of (3.3) on the stability boundary 1x1 = 1 is considerably more complicated. However, a number of results are available and they indicate that the attractor of the sequence {Y~,}~~~+ can be an exceedingly complicated object, as can be seen in Fig. 5. Unfortunately, virtually all (bar numerical experiments) that is known in precise mathematical terms with regard to the attractor refers to the special case X = - 1, which does not correspond to a realistic numerical method (since 9%= 1 + ha + 0(h2)). We expand each yn in powers of p, yn = C-V
+ ~a,,&, e=i
n E Z+.
304
A. lserles /Applied
Numerical Mathematics
24 (1997) 29.5-308
r
c
3 I L
X = -1, p
zz fi
X=i,p=$
Fig. 5. Attractors of the soldon sequence of (3.3) for different values of X and ,u, sketched in the plane (Re yn, Im yn) Let
US commence
y1 = -1
+2/i,
y2 = +1
-2p
by examining
the first few values
+2/_L2,
y3 =-I
++*,
y4 =+I
-2/L*
+2/x3,
-2/L*
+2p3,
y5 = -1
+2/i
35 = +I
-2/L
+2p3,
y7=-1
+G3
>
y/8 =+I
-2p3
+ 2y4,
y9 = -1
+2/l
-2p3 + 2p4,
y10 = +1
-2p
+2/L* -2p3 + 2p4 >
Y/II = -I
+2p2
-2p3 + 2p4,
Y12 = +I
-2/L*
-I- 2p?.
of pn,
A. Iserles /Applied Numerical Mathematics 24 (1997) 29.5-308
305
Obviously, we can discern a pattern! Firstly, it appears that an,e E (0, +2, -2) and, secondly, there is clearly a binary structure to the sequence. It is not difficult to show that this is indeed the case and that the explicit form of the coefficients is
0, “2e+Im+j.e
=
j =0,1,...,2e-l j=2e-l,...,2’j =2$...>3~2~-‘j=3x2”-‘....>2e+‘-l
2, -2, 0:
- 1: 1. 1.
for all e E N and m E Z+ [16]. Let A, be the w-set of the sequence {Yn}nEZ+. In other words, A, is the closure of all the accumulation points of the sequence. The explicit form of the coefficients allows to evaluate specific members of A,. For example, as long as jl_~J< 1 it is true that lim Y(CI)/X T.403
= -G
1-P
and
We present next an explicit evaluation 2’-1=2@l,+j
J&y2+1)/3
=
1-P
1 +pL’
of yzT_ 1. Since for all e = 1,2, . . . , T - 1 it is true that
wherem=2r-e-1-l,
j12~~l-1,
we deduce that ~~_l,e = 0 in that range. Moreover, for ! = T we have m = 0, j = 2’ - 1, whence CQ _ i,r = 2, and trivially Q_ 1,e = 0 for 1 3 T + 1. Therefore ~2~1 = - 1 + 2~’ and lim ~7-1
r+x
= -1.
It is possible to prove that, in general, each A, can be inscribed into a complex disc of radius (1 - r_ll/(l - 1~~1)[16]. It is known that A0 = {-l,+l} and A l/2 = [- 1, 11. Perhaps more surprising result is that every element of the sequence {Y%}~,=~+is in Ap. In order to prove this we pick an arbitrary n E Z+ and let rnk := n+2 Its-r, k E Z+. Obviously, the mk’s form a pairwise-distinct sequence and, provided that k is large enough, binary expansions of n and rnk differ only in the kth digit. It is easy to prove that arnlc,e = an,e, ! = 0, 1. . . ) k - 1, whence IYn-YlnkI s25p,~=3po B=k
for all Ir_~l< 1 [16]. As a matter of fact, unless CL= 0, the set A, is of an infinite cardinality. Moreover, it follows at once from the definition that there exists a countable number of components of A, in every nontrivial neighbourhood of every element of this set. Fig. 6 displays the sets A, for four choices of p E C. It is indicative of fractal sets, as is the above analysis. We proceed next to identify A,, p # 0, with a different mathematical construct, thereby proving that our suspicion is well-founded: the attractors indeed possess a fractional Hausdorff dimension. Letting X = -1 in (3.3), we eliminate pyn and pyn+i from these expressions. The first equation can be rewritten in the form ;PL(Y2n
PLY72 =
and substitution Y2nfl
-
+ Y2n+l)
in the second equation yields
Y2n =
y2n-1
-
Y2n-2.
306
A. Iserles /Applied
Numerical Mathematics
24 (1997) 295-308
I
0
0 p = +(l$i)
$i
,u =
1
1
0 -1
-10
0i? 4 $
I
J
0
p=
p=$-ai
&(l+
i)
Fig. 6. w-sets A, for different values of p E &1, displayed in the plane (Red,,
Consequently,
Imd,).
by induction,
Y2n+l - Y2n = Yl - Yo = 41
- p).
Adding and subtracting the first equation results in Y2n = PYn + (1 - P),
n E N,
Y2nSl
n E z+,
= PYn -
respectively. fl(4
(1 - P),
(3.5)
Letting := W + (1 - 4,
the recurrence
f2(4
:=
llz
-
(1
-
4,
z E c,
(3.5) becomes
Y2n = fl (Yn>,
Y2nfl
= f2(Yn).
This, in turn, has the same attractor as a dynamical system which is obtained when for each n E Z4 we obtain yn+t from yn by iterating either ft or f2 with equal probability [l]. The attractor .A, is thus identified with a probabilistic mixture of Julia sets, whose properties have been a subject to a thoroughgoing mathematical analysis in the last decade.
A. Iserles /Applied Numerical Mathematics 24 (1997) 295-308
307
To state that the present body of knowledge on the discretized pantograph equation is incomplete constitutes an understatement. Insofar as asymptotic stability is concerned, the method of proof extends just to Q which is a reciprocal of an integer and on the stability boundary we can explain only the nonrealistic case X = - 1. Yet, computer experimentation indicates that our results represent general patterns. It is fair to formulate the conjecture that, subject to further step-size restrictions in the neutral case, the discretization of (2.2) by an arbitrary multistep method is stable whenever (bl < la( and ha belongs to the linear stability domain of the underlying ODE scheme. Moreover, computer experiments seem to indicate that values of ha on the boundary of the linear stability domain lead to fractal attractors. Only limitations of our current methods of proof and of our understanding stand between us and a verification (or refutation) of these conjectures.
References [l] M.F. Barnsley and S.G. Den&o, Rational approximation of fractals, in: P.R. Graves-Morris, E.B. Saff and R.S. Varga, eds., Rational Approximation and Interpolation, Lecture Notes in Mathematics 1105 (Springer, Berlin, 1984) 73-88. [2] A. Bellen, N. Guglielmi and L. Torelli, Asymptotic stability properties of theta-methods for the pantograph equation, Quaderno Matematico 370 de1 Dipartimento di Scienze Matematiche dell’Universit8 di Trieste (1996). [3] H. Brunner, On the discretization of differential and Volterra integral equations with variable delay, Technical Report, Memorial University of Newfoundland (1996). [4] M.D. Buhmann and A. Iserles, On the dynamics of a discretized neutral equation, ZMA J. Numer. Anal. 12 (1992) 339-363. [5] M.D. Buhmann and A. Iserles, Stability of the discretized pantograph differential equation, Math. Comp. 60 (1993) 575-589. [6] M.D. Buhmann, A. Iserles and S.P. Norsett, Runge-Kutta methods for neutral differential equations, in: R.P. Agarwal, ed., Contributions in Numerical Mathematics, World Scientific Series in Applicable Analysis (World Scientific, Singapore, 1993) 85-98. [7] J. Carr and J. Dyson, The functional differential equation y’(x) = ay(Xz)+by(z), Proc. Roy. Sot. Edinburgh Sect. A 74 (1974/1975) 165-174. [8] J. Carr and J. Dyson, The matrix differential equation y’(x) = Ay(Xz) +By(z), Proc. Roy. Sot. Edinburgh Sect. A 75 (1975/1976) 5-22. of solutions of functional and differential-functional equations with several [91 G. Derfel, Behaviour transformations of the independent variable, Ukrai: Mat. Zh. 34 (1982) 350-356. [lo] G. Derfel, N. Dyn and D. Levin, Generalized refinement equation and subdivision processes, J. Approx. Theory 80 (1995) 100-111. [l l] G. Derfel and A. Iserles, The pantograph equation in the complex plane, DAMTP Technical Report 1996/NAOl (1996); also: J. Math. Anal. Appl., submitted. [12] L. Fox, D.F. Mayers, J.R. Ockendon and A.B. Tayler, On a functional-differential equation, J. Inst. Math. Appl. 8 (1971) 271-307. [13] K.J. in ‘t Hout, The stability of f&methods for systems of delay differential equations, Ann. Numer. Math. 1 (1994) 323-334. [ 141 A. Iserles, On the generalized pantograph functionaldifferential equation, European .I. Appl. Math. 4 (1993) l-38. [15] A. Iserles, On nonlinear delay differential equations, Trans. Amel: Math. Sot. 344 (1994) 441-477.
308
A. Iserles /Applied Numerical Mathematics 24 (1997) 295-308
[ 161 A. Iserles, The asymptotic behaviour of certain difference equations with proportional delay, Comput. Math. Appl. 28 (1994) 141-152. [17] A. Iserles and Y. Liu, On neutral functional-differential equations with proportional delays, DAMTP Technical Report 1993/NA3 (1993); also: J. Math. Anal. Appl., to appear. [18] A. Iserles and Y. Liu, On pantograph integro-differential equations, J. Integral Equations Appl. 6 (1994) 213-237. [19] T. Kato and J.B. McLeod, The functional-differential equation y’(z) = ay(Xx) + by(z), Bull. Amex Math. Sot. 77 (1971) 891-937. [20] Y. Liu, Stability of Q-methods for neutral functional-differential equations, Numel: Math. 70 (1995) 473485. [21] Y. Liu, On functional differential equations with proportional delays, Ph.D. Dissertation, University of
Cambridge ( 1996). [22] S. Skorik and V. Spiridonov, Self-similar potentials and the q-oscillator algebra at roots of unity, iAt. Math. Phys. 28 (1993) 59-74.