Found Com put M ath (2018) 18:1299-1332
https://doi.org/10.1007/s10208-017-9369-5
FOUNDATIONS o f
COMPUTATIONAL MATHEMATICS
The Journal of the Society for the Foundations of Computational Mathematics CrossMark
Algorithm for Rigorous Integration of Delay Differential Equations and the Computer-Assisted Proof of Periodic Orbits in the M ackey-Glass Equation
R o b e r t S z c z e lin a1 • P i o t r Z g lic z y iis k i2
Received: 4 July 2016 / Revised: 4 September 2017 / Accepted: 7 September 2017 / Published online: 30 October 2017
© The Author(s) 2017. This article is an open access publication
A b s t r a c t W e p re s e n t an a lg o rith m fo r th e rig o ro u s in te g ra tio n o f d e la y d iffe re n tia l e q u a tio n s (D D E s) o f th e fo rm x ' ( t ) = f ( x ( t — t), x ( t )). A s an a p p lic a tio n , w e g iv e a c o m p u te r-a s s is te d p r o o f o f th e e x is te n c e o f tw o a ttra c tin g p e rio d ic o rb its (b e fo re an d a fte r th e first p e rio d -d o u b lin g b ifu rc a tio n ) in th e M a c k e y - G la s s e q u a tio n .
K e y w o r d s C o m p u te r-a s s is te d p ro o fs ■ D e la y d iffe re n tia l e q u a tio n s ■ P e rio d ic o rb it ■ T o p o lo g ic a l m e th o d s ■ In te rv a l arith m e tic
M a t h e m a t i c s S u b je c t C la s s if ic a tio n 3 4 K 1 3 ■ 6 5 G 3 0 ■ 6 5 Q 2 0
1 Introduction
T h e g o a l o f th is p a p e r is to p re s e n t an a lg o rith m fo r th e rig o ro u s in te g ra tio n o f d e la y d iffe re n tia l e q u a tio n s (D D E s) o f th e fo rm
Communicated by Hans Munthe-Kaas.
B Robert Szczelina robert.szczelina@uj.edu.pl Piotr Zgliczynski zgliczyn@ii.uj.edu.pl
1 Malopolska Center of Biotechnology, Jagiellonian University, Gronostajowa 7a, 30-374 Krakow, Poland
2 Institute of Computer Science and Computational Mathematics, Jagiellonian University, Lojasiewicza 6, 30-348 Krakow, Poland
X ( t ) = f ( x ( t — t), x ( t ) ) , x e R (1)
<1 Springer u L ip j
w h e re 0 < t e R .
D e s p ite its a p p a re n t sim p licity , E q . ( 1) g e n e ra te s all k in d s o f p o s s ib le d y n a m ic a l b e h a v io rs : fro m sim p le sta tio n a ry so lu tio n s to c h a o tic a ttra c to rs. F o r e x a m p le , th is h a p p e n s fo r th e w e ll-k n o w n M a c k e y - G la s s e q u a tio n :
x (t — t) „
x ( t ) = p ■---y ■ x ( t ), x e R , (2) 1 + x n (t — t)
fo r w h ic h n u m e ric a l e x p e rim e n ts sh o w th e e x is te n c e o f a series o f p e rio d -d o u b lin g b ifu rc a tio n s w h ic h le a d to th e c re a tio n o f an a p p a re n t c h a o tic a ttra c to r [ 16, 17]. L a te r in th e p a p e r, w e w ill a p p ly o u r rig o ro u s in te g ra to r to th is e q u a tio n .
T h e re are m a n y im p o rta n t w o rk s th a t e sta b lish th e e x is te n c e an d th e s h a p e o f a (g lo b a l) a ttra c to r u n d e r v a rio u s a ss u m p tio n s on f in E q . ( 1) . M u c h is k n o w n a b o u t sy ste m s o f th e fo rm x = —^ x ( t ) + f ( x ( t — 1)) w h e n f is stric tly m o n o to n ic , e ith e r p o s itiv e o r n e g a tiv e [9 ]. L e t u s m e n tio n h e re a few d e v e lo p m e n ts in th is d ire c tio n . M a lle t-P a re t a n d S e ll u s e d d is c re te L y a p u n o v fu n c tio n a ls to p ro v e a P o in c a r e -B e n d ix s o n -ty p e th e o re m fo r sp e c ia l k in d o f m o n o to n e sy s te m s [ 19].
K ris z tin e t al. h a v e c o n d u c te d a th o ro u g h stu d y on sy ste m s h a v in g a m o n o to n e p o s itiv e fe e d b a c k , in c lu d in g stu d ie s on th e c o n d itio n s n e e d e d to o b ta in th e s h a p e o f a g lo b a l a ttra c to r; se e [11] an d re fe re n c e s th e re in . In th e c a s e o f a m o n o to n ic p o sitiv e fe e d b a c k f a n d u n d e r so m e a ss u m p tio n s on th e sta tio n a ry so lu tio n s, K ris z tin a n d Vas p ro v e d th a t th e re e x is t la rg e a m p litu d e slo w ly o s c illa to ry p e rio d ic so lu tio n s (L S O P s) w h ic h re v o lv e a ro u n d m o re th an o n e s ta tio n a ry so lu tio n . T o g e th e r w ith th e ir u n s ta b le m a n ifo ld s, c o n n e c tin g th e m w ith th e c la s sic a l sp in d le -lik e stru c tu re , th e y c o n s titu te th e fu ll g lo b a l a ttra c to r fo r th e sy s te m [ 10]. In a re c e n t w o rk , V as sh o w e d th a t f m ay b e c h o se n su ch th a t th e stru c tu re o f th e g lo b a l a ttra c to r m a y b e a rb itra rily c o m p lic a te d (c o n ta in in g an a rb itra ry n u m b e r o f u n s ta b le L S O P s ) [3 0 ].
L a n i-W a y d a a n d W a lth e r w e re a b le to c o n s tru c t sy ste m s o f th e fo rm x = f ( x ( t — 1)) fo r w h ic h th e y p ro v e d th e e x is te n c e o f tra n sv e rs a l h o m o c lin ic tr a je c tory, a n d a h y p e rb o lic se t on w h ic h th e d y n a m ic s are c h a o tic [ 13].
S rz e d n ic k i a n d L a n i-W a y d a p ro v e d , b y th e u s e o f th e g e n e ra liz e d L e fs h e tz fixed p o in t th e o re m , th e e x is te n c e o f m u ltip le p e rio d ic o rb its an d th e e x is te n c e o f c h a o s for so m e p e rio d ic , to o th -s h a p e d (p ie c e w ise lin e a r) f [ 12].
T h e re s u lts fro m [ 10, 12, 13, 3 0 ], w h ile im p re ss iv e , are e s ta b lis h e d fo r fu n c tio n s w h ic h are c lo s e to p ie c e w is e affin e o n es. T h e a u th o rs o f th e s e w o rk s c o n s tru c t e q u a tio n s w h e re an in te re stin g b e h a v io r a p p e a rs ; h o w ev er, it is n o t c le a r h o w to a p p ly th e ir te c h n iq u e s fo r so m e w e ll-k n o w n e q u a tio n s.
In re c e n t y e a rs , th e re a p p e a re d m a n y c o m p u te r-a s s is te d p ro o fs o f v a rio u s d y n a m ic a l p ro p e rtie s fo r o rd in a ry d iffe re n tia l e q u a tio n s an d (d issip a tiv e ) p a rtia l d iffe re n tia l e q u a tio n s b y an a p p lic a tio n o f a rg u m e n ts fro m th e g e o m e tric th e o ry o f d y n a m ic a l sy s te m s p lu s th e rig o ro u s in te g ra tio n ; see, fo r e x a m p le , [2 , 7 , 2 0 , 2 9 , 3 2 , 3 6 ] a n d r e f e re n c e s th e re in . B y th e c o m p u te r-a s s is te d p ro o f, w e u n d e rs ta n d a c o m p u te r p ro g ra m w h ic h rig o ro u s ly c h e c k s a ss u m p tio n s o f a b s tra c t th e o re m s. T h is p a p e r is an a tte m p t to e x te n d th is a p p ro a c h to th e c a s e o f D D E s b y c re a tin g a rig o ro u s fo rw a rd -in -tim e in te g ra tio n s c h e m e fo r E q . ( 1) . B y th e rig o ro u s in te g ra tio n w e u n d e rs ta n d a c o m p u te r p ro c e d u r e w h ic h p ro d u c e s rig o ro u s b o u n d s fo r th e tru e so lu tio n . In th e c a s e o f D D E s,
F0C71
1
Springeru
Found Com put M ath (2018) 18:1299-1332 1301
th e in te g ra to r sh o u ld re fle c t th e fa c t th at, a fte r th e in te g ra tio n tim e lo n g e r th a n the d e la y t , th e so lu tio n b e c o m e s sm o o th er, w h ic h g iv es th e c o m p a c tn e ss o f th e ev o lu tio n o p e ra to r. H a v in g an in te g ra to r, o n e sh o u ld b e a b le to d ire c tly a p p ly sta n d a rd to o ls fro m d y n a m ic s su c h as P o in c a re m a p s, v a rio u s fix ed p o in t th e o re m . In th is p ap er, as an a p p lic a tio n , w e p re s e n t c o m p u te r-a s s is te d p ro o fs o f th e e x is te n c e o f tw o stab le p e rio d ic o rb its fo r M a c k e y - G la s s e q u a tio n ; h o w ev er, w e d o n o t p ro v e th a t th e s e o rb its a re a ttra c tin g .
T h e re are se v e ra l p a p e rs th a t d e a l w ith c o m p u te r-a s s is te d p ro o fs o f p e rio d ic s o lu tio n s to D D E s [8 , 14, 3 3 ], b u t th e a p p ro a c h u s e d th e re is v e ry d iffe re n t fro m o u r m eth o d . T h e s e w o rk s tra n s fo rm th e q u e stio n o f th e e x is te n c e o f p e rio d ic o rb its in to a b o u n d a ry v a lu e p ro b le m (B V P ), w h ic h is th e n so lv e d b y u s in g th e N e w to n -K a n to ro v ic h th e o re m [8 , 14] or th e lo c a l B ro u w e r d e g re e [3 3 ]. It is c le a r th a t th e rig o ro u s in te g ra tio n m ay b e u s e d to o b ta in m o re d iv e rse sp e c tru m o f re s u lts . T h e re are a lso sev eral in te re stin g re s u lts th a t a p p ly rig o ro u s n u m e ric a l c o m p u ta tio n s to so lv e p ro b le m s fo r D D E s [3 , 4 ], b u t th ey d o n o t re ly on th e rig o ro u s, fo rw a rd -in -tim e in te g ra tio n o f D D E s.
T h e re s t o f th e p a p e r is o rg a n iz e d as fo llo w s. S e c tio n 2 d e sc rib e s th e th e o ry and a lg o rith m s fo r th e in te g ra tio n o f E q . ( 1) . S e c tio n 3 d efin es th e n o tio n o f th e P o in c a re m a p a n d d is c u s s e s c o m p u ta tio n o f th e P o in c a re m ap u s in g th e rig o ro u s in teg rato r.
S e c tio n 4 p re s e n ts an a p p lic a tio n o f th e m e th o d to p ro v e th e e x is te n c e o f tw o sta b le p e rio d ic o rb its in th e M a c k e y - G la s s e q u a tio n (E q. (2 ) ). H e re , w e in v e s tig a te c a s e for n = 6 (b e fo re th e first p e rio d -d o u b lin g b ifu rc a tio n ) a n d fo r n = 8 (a fte r th e first p e rio d -d o u b lin g b ifu rc a tio n ). To th e b e s t o f o u r k n o w le d g e , th e s e are th e first rig o ro u s p ro o fs o f th e e x is te n c e o f th e s e o rb its. T h e p re s e n te d m e th o d h a s b e e n a lso su c c e s sfu lly u s e d b y th e first a u th o r to p ro v e th e e x is te n c e o f m u ltip le p e rio d ic o rb its in so m e o th er n o n lin e a r D D E s [2 5 ].
1.1 N o ta tio n
W e u s e th e fo llo w in g n o ta tio n . F o r a fu n c tio n f : R ^ R , b y f (k), w e d e n o te the k t h d e riv a tiv e o f f . B y f [k], w e d e n o te th e te rm 1 ■ f (k). In th e c o n te x t o f p ie c e w ise sm o o th m a p s b y f (k)( t - ) an d f (k)( t + ) , w e d e n o te th e o n e -s id e d d e riv ativ es f w .r.t.
t .
F o r F : R m ^ R n b y D F ( z ) , w e d e n o te th e m a trix ( ( z ) ) .
\ dxj ) i e [1,...,n}, j
F o r a g iv en se t A , b y cl (A ) an d in t (A ), w e d e n o te th e c lo s u re a n d in te rio r o f A, re s p e c tiv e ly (in a g iv en to p o lo g y , e.g ., d e fin e d b y th e n o rm in th e c o n s id e re d B a n a c h sp ace).
L e t A = n n = l [ai, b i ] fo r a i < b i , a i , bi e R . T h e n , w e c a ll A a n in t e r v a l set (a p ro d u c t o f c lo s e d in te rv a ls in R n ). F o r a n y A c R n , w e d e n o te b y h u l l ( A ) a m in im a l in te rv a l set, su ch th a t A c h u l l ( A ) . I f A c R is b o u n d e d th e n h u l l ( A ) = [in f(A ), su p (A )]. F o r sets A c R , B c R , a e R an d fo r so m e b in a ry o p e ra tio n o : R x R ^ R w e d e f in e A o B = [a o b : a e A , b e B } a n d a o A = A o a = [a } « A . A n a lo g o u sly , fo r g : R ^ R an d a se t A e R w e d e fin e g ( A ) = { g ( a ) | a e A}.
F o r v e R n b y n i v fo r i e [1 , 2 , . . . , n } , w e d e n o te th e p ro je c tio n o f v o n to the i th c o o rd in a te . F o r v e c to rs u , v e R n b y u ■ v, w e d e n o te th e sta n d a rd sc a la r p ro d u ct:
u ■ v = J 2 n= 1 n i v ■ n i u
FoCTI
1
Springer uW e d e n o te b y C r (D , R ) th e sp a c e o f all fu n c tio n s o f c la s s C r o v er a c o m p a c t set D c R , e q u ip p e d w ith th e s u p re m u m C r n o rm : ||g || = Y ,r = 0 su p x eD |g (l\ x ) |. In c a s e D = [—t, 0 ], w h e n t is k n o w n fro m th e c o n te x t, w e w ill w rite C k in s te a d o f C k ( [ -t, 0 ], R ).
F o r a g iv en fu n c tio n x : [ —1, a ) ^ R , a e R + U ( c ^ j fo r an y t e [0, a ) w e d e n o te b y x t a fu n c tio n su ch th a t x t (s ) = x ( t + s ) fo r all s e [- 1 , 0].
W e w ill o ften u s e a sy m b o l in s q u a re b ra c k e ts, e .g ., [r], to d e n o te a se t in R m . U s u a lly it w ill h a p p e n in fo rm u la s u s e d in a lg o rith m s, w h e n w e w o u ld lik e to stress th e fa c t th a t a g iv en v a ria b le re p re s e n ts a set. I f b o th v a ria b le s r a n d [r] are u se d s im u lta n e o u sly , th e n u s u a lly r re p re s e n ts a v a lu e in [r]; h o w ev er, th is is n o t im p lie d b y d e fa u lt an d it w ill b e alw ay s sta te d e x p lic itly . P le a s e n o te th a t th e n o ta tio n [r] does n o t im p o s e th a t th e set [r] is o f an y p a rtic u la r sh ap e, e .g ., an in te rv a l b o x . W e w ill a lw a y s e x p lic itly sta te i f th e se t is an in te rv a l b o x .
F o r an y se t X b y m i d ( X ), w e d e n o te th e m id p o in t o f h u l l ( X ) a n d b y d i a m ( X ) the d ia m e te r o f h u l l ( X ).
1 .2 B a s ic P r o p e r t i e s o f S o lu tio n s to D D E s
F o r th e c o n v e n ie n c e o f th e re a d e r, w e re c a ll (w ith o u t p ro o f s ) sev eral c la s sic a l re s u lts fo r D D E s [5 ].
W e d e fin e th e s e m if lo w y a ss o c ia te d w ith E q . ( 1) by:
y : R + x C 0 ( [ —t, 0 ], R ) 9 ( t , f ) ^ x f e C 0 ( [ —t, 0 ], R ) . (3)
w h e re x f : [—t, a f ) ^ R is a so lu tio n to a C a u c h y p ro b le m :
(4) I x = f ( x ( t - t), x ( t ) ) , t > 0,
| x ( t ) = f ( t ), t e [ —t , 0],
fo r a m a x im a l a f e R + U ( c » j su ch th a t th e so lu tio n e x is ts fo r all t < a f .
L e m m a 1 [C o n tin u o u s (lo c a l) sem iflo w ] I f f is (locally) L ipshitz, th e n y is a (local) c o n ti n u o u s s e m if lo w o n C0( [ —t, 0 ], R ).
L e m m a 2 [S m o o th in g p ro p e rty ] A s s u m e f is o f c la s s C m, m > 0. L e t n e N b e g iv e n a n d l e t t > n ■ t. I f x 0 e C0 t h e n x t = y ( t , x 0) is o f c la s s a t le a s t C min(m+1,n).
T h e sm o o th in g o f so lu tio n s g iv es ris e to so m e in te re stin g o b je c ts in D D E s [3 1 ].
A s s u m e fo r a w h ile th a t f e C “ . T h e n fo r an y n > 0, th e re e x is ts a se t (in fa c t a m a n ifo ld ) M n c C n , su c h th a t M n is fo rw a rd in v a ria n t u n d e r y .
It is e a s y to se e th a t fo r n = 1 w e have:
M 1 : = { x e C1 | x '(0—) = f ( x ( —t), x (0) ) } ,
a n d th e c o n d itio n s fo r M n w ith n > 1 c a n b e sim p ly o b ta in e d b y d iffe re n tia tin g b o th sid es o f ( 1) . W e fo llo w [3 1 ] a n d w e c a ll M n a C n s o l u ti o n m an ifo ld .
FoCJI
1
Springeru
Found Com put M ath (2018) 18:1299-1332 1303
N o tic e th a t M n c M k fo r k < n a n d ^ ( Kt, ■) : M n ^ M n+ k.
2 Rigorous Integration of DDEs
T h is se c tio n is a re o rg a n iz e d e x c e rp t fro m th e P h .D . d is se rta tio n o f th e first au th o r (R o b e rt S z c z e lin a ). A d e ta ile d a n a ly sis o f re s u lts fro m n u m e ric a l e x p e rim e n ts w ith th e p ro p o s e d m e th o d s, m o re e la b o ra te d e s c rip tio n o f th e a lg o rith m s, a n d d e ta ile d p s e u d o -c o d e s o f th e ro u tin e s c a n b e fo u n d in th e o rig in a l d is s e rta tio n [2 4 ].
2 .1 F in i te R e p r e s e n t a ti o n o f “ S u ff ic ie n tly S m o o th ” F u n c ti o n s
H e re , w e w o u ld lik e to p re s e n t th e b a sic b lo c k s u s e d in th e a lg o rith m fo r th e rig o ro u s in te g ra tio n o f E q . ( 1) . T h e id e a is to im p le m e n t th e T ay lo r m e th o d fo r E q . ( 1) b a se d on th e p ie c e w is e p o ly n o m ia l re p re s e n ta tio n o f th e so lu tio n s p lu s a re m a in d e r term . W e w ill w o rk on th e e q u a lly sp a c e d g rid an d w e w ill fix th e step siz e o f th e T ay lo r m e th o d to m a tc h th e s e le c te d g rid .
R e m a r k 1 In th is se c tio n , fo r th e sa k e o f s im p lic ity o f p re s e n ta tio n , w e a ss u m e th at
t = 1. A ll c o m p u ta tio n s c a n b e e a sily re d o n e fo r an y d e la y t .
W e also a ss u m e th a t r.h .s. f o f E q . ( 1) is “ su ffic ie n tly sm o o th ” fo r v a rio u s e x p re s sio n s to m a k e sen se. T h e class o f f in ( 1) re s tric ts th e p o s s ib le o rd e r o f th e T ay lo r m e th o d th a t c a n b e u s e d in o u r a lg o rith m s, th a t is, i f f is o f class C n , th e n w e can u se T a y lo r m e th o d o f o rd e r a t m o s t n . T h e re fo re , th o ro u g h th e p a p e r it c a n b e a ss u m e d th a t f e C “ . T h is is a re a s o n a b le a s s u m p tio n in th e c a se o f a p p lic a tio n s o f c o m p u te r- a ss is te d p ro o fs w h e re r.h .s. o f e q u a tio n s are u s u a lly p re s e n te d as a c o m p o s itio n o f e le m e n ta ry fu n c tio n s. T h e M a c k e y - G la s s e q u a tio n (2 ) is a g o o d e x a m p le (a w a y fro m x = - 1).
W e fix tw o in te g e rs n > 0 a n d p > 0 a n d w e se t h = p .
D e f in itio n 1 B y C np , w e d e n o te th e se t o f all fu n c tio n s g : [ —1, 0 ] ^ R su ch th at, fo r i e {1, . . . , p } , w e have:
• g is (n + 1) -tim es d iffe re n tia b le on ( —i ■ h , —i ■ h + h ) ,
• g (k) ( —i ■ h + ) e x is ts fo r all k e { 0 , . . . , n + 1} a n d lim g (k)( —i ■ h + %) =
£ ^0+ g (k)( —i ■ h + ) ,
• g^(n+ 1') is c o n tin u o u s an d b o u n d e d on ( —i ■ h , —i ■ h + h ) .
F ro m n o w on, w e w ill a b u se th e n o ta tio n an d w e w ill d e n o te th e rig h t d eriv ativ e g (k)( —i ■ h + ) b y g (k)( —i ■ h ) u n le s s e x p lic itly sta te d o th e rw ise . T h e sa m e h o ld s for g [k]( —i ■ h + ) . U n d e r th is n o ta tio n , it is c le a r th a t w e ca n re p re s e n t g e C np b y a p i e c e w i s e T aylor e x p a n sio n on e a c h in te rv a l [—i ■ h , —i ■ h + h ) . F o r t = —i ■ h + £
a n d 1 < i < p , 0 < £ < h w e c a n w rite:
n
g ( t ) =
J 2
g [k]( —i ■ h ) ■ £k + g [n+ 1] ( —i ■ h + £ ( £ ) ) ■ £n+ 1, (5) k=0FoCTI u
< 1 S p rin g e r
w ith %(e) e [0, h ].
In o u r a p p ro a c h , w e sto re th e p ie c e w ise T a y lo r e x p a n sio n as a fin ite c o lle c tio n o f c o e ffic ie n ts g [k]( - i ■ h ) an d in te rv a l b o u n d s on g [n+1](-) o v er th e w h o le in terv al [—i ■ h , - i ■ h + h ] fo r i e { 1 , . . . , p } . O u r a lg o rith m fo r th e rig o ro u s in te g ra tio n o f ( 1) w ill th e n p ro d u c e rig o ro u s b o u n d s on th e so lu tio n s to (1) fo r in itia l fu n c tio n s d efin ed b y su c h p ie c e w is e T ay lo r e x p a n sio n .
P le a s e n o te th a t w e are u s in g h e re a w o rd f u n c t i o n s in s te a d o f a sin g le f u n c t i o n , as, b e c a u s e o f th e b o u n d s on g [n+ 1] ov er in te rv a ls [—i ■ h , —i ■ h + h], th e fin ite p ie c e w ise T a y lo r e x p a n sio n d e sc rib e s an in fin ite se t o f fu n c tio n s in g e n e ra l. T h is m o tiv a te s the fo llo w in g d e fin itio n s.
D e f in itio n 2 L e t g e C p an d le t I : { 1 , . . . , p } x { 0 , . . . , n } ^ { 1 , . . . , p ■ (n + 1)}
b e an y b ije c tio n .
A m i n i m a l (p , n ) - r e p r e s e n t a ti o n o f g is a p a ir g = (a, B ) su c h th a t
a e '■(«+i)+ib c R p is an in te rv a l set
n \ a = g ( 0 ) n i B =
K1+ m ,k)a = g [k]( —i h + ) for 1 < i < p , 0 < k < n
in f g [n+ 1] ( —i h + %), su p g [n+ 1] ( —i h + %)
%e(0,h) %e(0,h)
P le a s e n o te , th a t th e in d e x fu n c tio n I sh o u ld b e sim p ly u n d e rs to o d as an o rd e rin g in w h ic h , d u rin g c o m p u ta tio n s , w e sto re c o e ffic ie n ts in a fin ite -d im e n sio n a l v e c to r— its p re c is e d e fin itio n is o n ly im p o rta n t fro m th e p ro g ra m m in g p o in t o f view , se e S ect. 2.4 fo r a p a rtic u la r e x a m p le o f I . S o, in th is p a p e r fo r th e o re tic a l c o n s id e ra tio n s, w e w o u ld lik e to u s e th e fo llo w in g g n o ta tio n in stead :
• g°,[°] : = n1a ,
• g i,[k] := n 1+i(i,k)a,
• gi,[n+1] := n B
W e c a ll g 1' ^ th e (i,k)th c oeffic ient o f th e r e p r e se n ta tio n an d g I,[” + 1] th e ith r e m a in d e r o f th e repr esen ta tio n . T h e in te rv a l set B is c a lle d th e r e m a i n d e r o f th e representa tio n.
W e w ill c a ll th e c o n sta n t m = p ■ (n + 2) + 1 th e size o f th e (p , n ) - re p re se n ta tio n . W h e n p a ra m e te rs n an d p are k n o w n fro m th e c o n te x t, w e w ill o m it th e m an d w e w ill c a ll g th e m i n i m a l r e p r e se n ta tio n o f g .
D e f in itio n 3 W e sa y th a t G c C np is a (p , n ) - f - s e t (o r (p , n )- fu n c tio n s set) i f th e re e x ists b o u n d e d se t [g] = (A , C ) c Rp ' (n+ 1 )+ 1 x R p = R m su ch th a t
G = | f e Cp | f c [ g ] fo r th e m in im a l (p , n ) -re p re se n ta tio n f o f f J .
A s th e se t [g] c o n ta in s th e m in im a l re p re s e n ta tio n o f f fo r an y f e G , w e w ill also say th a t [g] is a ( p , n ) - re p re s e n ta tio n o f G . W e w ill also u s e G an d [g] in te rc h a n g e a b ly a n d w e w ill w rite f e [g] fo r sh o rt, i f th e c o n te x t is clear.
F
0
C71
u<1 Springer
Found Com put M ath (2018) 18:1299-1332 1305
P le a s e n o te th a t th e m in im a l (p , n ) -re p re se n ta tio n g o f g d e fin e s (p , n ) - f - s e t G c C np , w h ic h , in g e n e ra l, c o n ta in s m o re th a n th e so le fu n c tio n g . A ls o , in g e n e ra l, fo r any ( p , n ) - f - s e t G th e re a re fu n c tio n s g e G w h ic h are d is c o n tin u o u s a t g rid p o in ts —i ■ h (s e e (5 )). S o m e tim e s, w e w ill n e e d to a ss u m e h ig h e r re g u la rity , th e re fo re w e define:
D e f in itio n 4 L e t G = [g] b e a (p , n )- f-s e t. T he C k - s u p p o r t o f G is d e fin e d as:
S u p p (k)([g]) := S u p p (k)( G ) : = G n C k .
F o r c o n v e n ie n c e w e also set:
S u p p ( [ g ] ) := S u p p ( G ) : = G.
P le a s e n o te th at, in g e n e ra l, C np d S u p p ( G ) = S u p p (0)( G ) c C0. It m a y also h a p p e n th a t S u p p (k)( G ) = 0 fo r n o n -e m p ty G even fo r k = 0.
N o w w e p re s e n t th re e sim p le facts a b o u t th e co n v e x ity o f th e su p p o rt sets. T h e se p ro p e rtie s w ill b e im p o rta n t in th e c o n te x t o f th e c o m p u te r-a ssiste d p ro o fs an d in an a p p lic a tio n o f T h e o re m 11 to ( p , n ) -f-se ts in S ect. 4 .
L e m m a 3 F or (p , n ) - r e p r e s e n t a ti o n s [j>], [ f ] c R m, m = p ( n + 2 ) + 1 , th e f o l l o w i n g s t a t e m e n t s h o l d tr u e :
- I f [g] c [ f ] , t h e n S u p p ( [ g ] ) c S u p p ( [ f ] ) .
- I f [g] is a c o n v e x s e t in R m, th e n S u p p ( [ g ] ) is a c o n v e x se t in C np .
- I f [g ] is a c o n v e x s e t in R m, th e n S u p p ( [ g ] ) n C k is a c o n v e x s e t f o r a n y k > 0.
W e o m it th e e a sy p ro o f.
To e x tra c t in fo rm a tio n on g [k] ^ — p + e ^ fo r an y i a n d k h a v in g o n ly in fo rm a tio n s to re d in a (p , n ) -re p re se n ta tio n , w e in tro d u c e th e fo llo w in g d efin itio n .
D e f in itio n 5 L e t ( p , n )- re p re s e n ta tio n g b e g iv en . W e d efin e
c f k](e) = W k ) ■ e l—k ■ g i,[l\
l=k ' '
fo r 0 < e < 1 , 1 < i < p a n d 0 < k < n + 1.
W e w ill o m it su b s c rip t g in c g [k] (e ) i f it is c le a r fro m th e co n te x t. T h e fo llo w in g le m m a fo llo w s im m e d ia te ly fro m th e T a y lo r fo rm u la , so w e sk ip th e p ro o f:
L e m m a 4 A s s u m e g e C np a n d its (p , n ) - r e p r e s e n t a ti o n g are giv en. T h e n f o r 0 <
e < p , 1 < i < p a n d 0 < k < n + 1
g M L + e ^ e c l' [k (e)
holds.
FoCTI u
<1 Springer
B e fo re p ro c e e d in g to th e p re s e n ta tio n o f th e in te g ra tio n p ro c e d u re , w e w o u ld lik e to d is c u ss th e p ro b le m o f o b ta in in g T a y lo r c o e ffic ie n ts o f a so lu tio n x to E q . ( 1) at a g iv en tim e t (w h e n e v e r th e y e x ist). F ro m E q . ( 1), w e h a v e (w e re m in d th at, at g rid p o in ts , b y th e d e riv a tiv e w e m e a n th e r ig h t d eriv ativ e):
d k —1
x (k)( t ) = d t k - 1 f ( x ( t — 1) ’ x ( t ))■
F o r e x a m p le , in c a s e o f k = 1, w e o b v io u s ly have:
x (1)( t ) = f ( x ( t — 1), x ( t ) ) ,
a n d in c a s e k = 2, b y a p p ly in g th e c h a in ru le , w e get:
x {2)( t ) = d - ( x ( t — 1), x ( t )) ■ x (1)(— 1) d z1
+ d - ( x ( t — 1), x ( t )) ■ x (1)(0) d z2
I f w e d e fin e a fu n c tio n F (i) : R4 ^ R as
d f d f
F (1) ( z 1 , z 2 ’ z 3 ’ z4) = — ( z 1 , z3) ■ z 2 + — ( z 1 , z3) ■ z4,
dz1 dz3
th e n w e se e th a t
x (2>( t ) = F(t) ( x ( t — 1), x (1)(t — 1), x ( t ), x (1) (t) ^ .
N ow , b y a re c u rs iv e a p p lic a tio n o f th e c h a in ru le , w e ca n o b ta in a fa m ily o f fu n c tio n s F f = {F(k) : R 2'(k+ 1) ^ R }keN su c h that:
x (k+ 1)( t ) = F (k) ( x ( — 1 + t ) , . . . ’ x (k)( — 1 + 1), x ( t ) , . . . ’ x (k)( t ) ) .
B y settin g
F [ k ] ( z 1 ,. . . ’ z2(k+1)) = k j F(k) (0! ■ z 1 , . . . ’ k !■ zk+1, 0! ■ zk+2 , . . . , k !■ z2-(k+1)), (6) w e c a n w rite sim ila r id e n tity in te rm s o f th e T a y lo r c o e ffic ie n ts x [k]:
x [k+ 1]( t ) = — ^ ■ F[k] ( x [0] ( —1 + 1) , . . . ’ x [k]( —1 + t ), x [0]( t ) , . . . ’ x [k]( t ) ) . (7)
k + 1 ^ /
A s w e are u sin g th e T ay lo r c o e ffic ie n ts in s te a d o f d e riv ativ es to re p re s e n t o u r ( p , n )- f-sets, th is n o ta tio n w o u ld b e m o re su ita b le to d e s c rib e c o m p u te r a lg o rith m s. F ro m n o w on, w e w ill also slig h tly a b u se th e n o ta tio n an d w e w ill d e n o te F[k] b y F [k]
FoCTI
<1 Springer
u
Found Com put M ath (2018) 18:1299-1332 1307
a n d F(k ) b y F (k). T h is is re a s o n a b le , sin ce, fo r a fu n c tio n F : R ^ R d e fin e d b y F ( t ) : = f ( x ( t — 1), x ( t )), w e have:
F (k)( t ) = F(k) ( x (0)( —1 + t ) , . . . , x (k)( —1 + 1), x (0) ( t ) , . . . , x (k)]( t ) ) , a n d
F [k]( t ) = F[k] ( x [0]( —1 + t ) , . . . , x [k]( —1 + t ), x [0] ( t ) , . . . , x [k\ t ) ) .
R e m a r k 2 T h e ta s k o f o b ta in in g fa m ily F f b y d ire c tly a n d a n a ly tic a lly a p p ly in g the c h a in ru le m a y s e e m q u ite te d io u s, esp e c ia lly , i f o n e w ill b e re q u ire d to su p p ly th is fa m ily as im p le m e n ta tio n s o f c o m p u te r p ro c e d u re s . It tu rn s o ut, th a t th is is n o t th e case fo r a w id e cla s s o f fu n c tio n s. In fact, o n ly th e r.h .s. o f E q . ( 1) n e e d s to b e im p le m e n te d a n d th e d eriv ativ es m a y b e o b ta in e d b y th e m e a n s o f th e au to m a tic d iffe re n tia tio n (A D ) [2 1 , 2 6 ]. W e u s e T ay lo r c o e ffic ie n ts x [k] to fo llo w th e n o ta tio n an d im p le m e n ta tio n o f A D in th e C A P D lib ra ry [ 1] w h ic h p ro v id e a se t o f rig o ro u s in te rv a l a rith m e tic ro u tin e s u s e d in o u r p ro g ra m s.
2 .2 O n e S te p o f th e I n t e g r a t i o n w ith F ix e d - S iz e S te p h = p
W e a re g iv e n (p , n ) - f - s e t x0 an d th e ta s k is to o b ta in x h— a ( p , n ) - f - s e t su c h th a t y ( h , x0) C x.h . W e w ill d e n o te th e p ro c e d u re o f c o m p u tin g x h b y I h , th a t is:
xh = Ih ( x0 ).
F irs t o f all, w e c o n s id e r h o w x0 a n d x h = I h (x 0) r e la te to e a c h other. T h e ir m u tu al a lig n m e n t is sh o w n in F ig . 1.
W e see th a t x ^ o v erlap w ith x0—1,[k], so th e y ca n b e sim p ly s h i f t e d to th e n ew re p re s e n ta tio n — w e c a ll th is p ro c e d u r e th e shi f t p a r t . O th e r c o e ffic ie n ts n e e d to b e e s tim a te d u s in g th e d y n a m ic s g e n e ra te d b y E q . ( 1) . W e c a ll th is p ro c e d u re th e f o r w a r d p a r t . T h is p ro c e d u re w ill b e d iv id e d in to th re e su b ro u tin es:
1. c o m p u tin g c o e ffic ie n ts x l ’[k fo r k e {1, . . . , n},
2. c o m p u tin g th e re m a in d e r x l ’[n + 1 ,
3. c o m p u tin g th e e s tim a te fo r x h (0) (s to re d in x ° ,[0]).
F o r w a r d P a r t : S u b r o u t i n e 1
T h is p ro c e d u re is im m e d ia te ly o b ta in e d b y a re c u rs iv e a p p lic a tio n o f E q . (7 ) :
D ,[0] _ =0,[0] ' h
= 1,[k]
= x,0
= - ■ f [k—1] x p,[0]
x' 0 , z.p,[k~1] x 1,[0]
, x' h , x lh [k—l ] \
w h e re 1 < k < n .
FoCTI u
<1 Springer
h 0
Fig. 1 A graphical presentation of the integrator scheme. We set n = 2 and p = 4.A (p, n) -representation is depicted as dots at grid points and rectangles stretching on the whole intervals between consecutive grid points. The dot is used to stress the fact that the corresponding coefficient represents the value at a given grid point. Rectangles are used to stress the fact that remainders are bounds for derivative over whole intervals.
Below the time line we have an initial (p, n)-representation. Above the time line, we see a representation after one step of size h = 1 . Black solid dots and gray rectangles represent the values we do not need to compute—this is the shift part. The forward part, i.e., the elements to be computed, are presented as empty dots and an empty rectangle. The doubly bordered dot represents the exact value of the solution at the time t = h = p (in practical rigorous computations it is an interval bound on the value). The doubly bordered empty rectangle is an enclosure for the n + 1th derivative on the interval [0, h]
F o C J l U
<1 Springer
Found Com put M ath (2018) 18:1299-1332 1309
F o r w a r d P a r t : S u b r o u t in e 2
T h is s u b ro u tin e c a n b e d e riv e d fro m th e m e a n v a lu e th e o re m . W e h a v e fo r e < h:
1 x (n +1) (e ) = ---1--- x (n+1) (0)--- + -1--- x (n+2) ( £) ■ e =
(n + 1)! W (n + 1)! W (n + 1)!
F [n] ( x ( —1) , . . . , x [n\ — 1), x ( 0 ) , . . . , x [n](0), ) (n + 1)
+ F [n+ 1] ( x ( —1 + £ ) , . . . , x [n+ 1] ( —1 + £), x ( £ ) , . . . , x [n+ 1]( £ ) ) ■ e
fo r so m e 0 < £ < e . L e t u s lo o k a t th e tw o te rm s th a t a p p e a r on th e r.h .s. o f E q . (8) . T h e q u e stio n is: C a n w e e s tim a te th e m b y h a v in g o n ly x0 an d a lre a d y c o m p u te d x ^ [k]
fro m S u b ro u tin e 1 ? L e t u s d is c u ss e a c h o f th e s e te rm s sep arately . B y L e m m a 4 , w e h a v e fo r 0 < k < n + 1:
x [k]( —1 + £) e c f0[k] ([0, h ])
M o re o v e r, b y D e fin itio n 2, w e k n o w th at:
x [k]( —1) e x p ,[k]
x [k](0) e x h1,[k]
S o th o s e te rm s c a n b e e a s ily o b ta in e d . T h e p ro b le m a p p e a rs w h e n it c o m e s to x (£) a n d x [k](£ ), fo r £ e (0, h )
A s s u m e fo r a m o m e n t th a t w e h a v e so m e a p rio ri e stim a te s fo r x ([0, h ]), i.e ., a se t Z c R su c h th a t x ([0, h ] ) C Z . W e c a ll th is se t th e ro u g h e n c lo s u re o f x on the in te rv a l [0, h ] . H a v in g ro u g h e n c lo s u re Z , w e c o u ld a p p ly E q . (7 ) (as in th e c a s e o f S u b ro u tin e 1) to o b ta in th e e stim a te s on x [k] ([0, h ] ) fo r k > 0. S o th e q u e stio n is:
h o w to find a c a n d id a te Z an d p ro v e th a t x ([0 , h ]) C Z ? T h e fo llo w in g le m m a gives a p ro c e d u r e to te s t th e later.
L e m m a 5 L e t Y c R b e a c lo s e d in te r v a l a n d le t x 0 b e a fu n c t io n d e fin e d o n [ —1, 0], A s s u m e th a t th e fo llo w in g h o ld s tr u e :
Z : = x0(0) + [0, h ] ■ f ( x0 ( [ —1, —1 + h ] ) , Y ) C i n t ( Y ). (8)
T h e n th e s o lu tio n x ( t ) o f E q, ( 1) w ith th e in itia l c o n d itio n x 0 e x is ts o n th e in te rv a l
[ 0, h ] a n d
x ([0, h ] ) C Z .
P r o o f W e c a n tre a t E q . ( 1) on th e in te rv a l [0, h ] as a n o n -a u to n o m o u s O D E o f the form :
x ' = f ( a ( t ), x ),
1
Springerw h e re a ( t ) = x ( t ) fo r t e [—1, —1 + h] is a k n o w n fu n c tio n . N o w th e c o n c lu s io n fo llo w s fro m th e p r o o f o f th e a n a lo g o u s th e o re m fo r O D E s. T h e p r o o f c a n b e fo u n d
in [3 5 ]. □
U sin g L e m m a 5 , a h e u ris tic ite ra tiv e a lg o rith m m a y b e d e s ig n e d su c h th a t it starts b y g u e s s in g an in itia l Y , a n d th en it ap p lie s E q . (8) to o b ta in Z . In a c a s e o f fa ilu re o f th e in c lu sio n , i.e ., Z c Y , a b ig g e r Y is ta k e n in th e n e x t ite ra tio n . P le a s e n o te th a t th is ite ra tio n m a y n ev er sto p o r p ro d u c e u n a c c e p ta b ly b ig Y , e s p e c ia lly w h e n th e ste p siz e h is larg e. F in d in g a ro u g h e n c lo s u re is th e o n ly p la c e in th e a lg o rith m o f th e in te g ra to r th a t c a n in fa c t fail to p ro d u c e a n y e stim a te s. In su ch a c ase, w e a re n o t ab le to p ro c e e d w ith th e in te g ra tio n , an d w e s ig n a liz e an error.
N o w w e c a n s u m m a riz e th e a lg o rith m fo r S u b ro u tin e 2 as fo llo w s:
c [k] : = c p ,[k]([0, h ] ) , k e ( 0 , . . . , n + 1} d [0] : = Z as in E q .(8) ,
d [k] : = 1 ■ F [k—1] ( c [0], . . . , c [k—1], d [0], . . . , d [k—1]) , k e (1, a * ____ 1 F [n] ( x p,[°] x pAn} -1,[0] -1, [n]) u := (n + 1) • c ^ 0 , . . . , A 0 , Ah , . . . , Ah j b* : = F [n+1] ( c [0] , . . . , c [n+1] , d [0], . . . , d [n+1]),
x,1,[n+1] : = a* + b* ■ [0, h ]
, n + 1}
R e m a r k 3 P le a s e n o te th a t th e te rm a* is c o m p u te d th e sa m e w a y as o th e r co efficien ts in S u b ro u tin e 1 an d th e ro u g h e n c lo s u re d o n o t in flu e n c e th is term . In fa c t th is is th e n + 1th d e riv a tiv e o f th e flow w .r.t. tim e . It is p o s s ib le to k e e p tra c k o f th o se c o e ffic ie n ts d u rin g th e in te g ra tio n an d a fte r p step s (fu ll d e la y ) th o s e c o e ffic ie n ts m a y b e u s e d to b u ild a (p , n + 1)-re p re s e n ta tio n o f th e s o lu tio n s— th is is a d ire c t re fle c tio n o f c o n se q u e n c e s o f L e m m a 2 .
T h is fa c t is also im p o rta n t fo r th e co m p a c tn e ss o f th e e v o lu tio n o p e ra to r— an e s s e n tia l p ro p e rty th a t allo w s fo r an a p p lic a tio n o f th e to p o lo g ic a l fixed p o in t th e o re m s in in fin ite -d im e n s io n a l sp aces.
F o r w a r d P a r t : S u b r o u t i n e 3
T h e la s t s u b ro u tin e o f th e fo rw a rd p a rt ca n b e sim p ly o b ta in e d b y u sin g D e fin itio n 2 a n d E q . (5 ) :
n
x h° ,[0] : = £ x1' ^ h k + x f :n+1] - h n+ 1. k =0
N o tic e th a t th e p o s s ib le in flu e n c e o f th e u s u a lly o v e re s tim a te d ro u g h e n c lo s u re Z is p re s e n t o n ly in th e la s t te rm o f th e o rd e r h n+ 1 , so fo r sm a ll h (la rg e e n o u g h p ), it sh o u ld n o t b e a p ro b le m .
T h e I n t e g r a t o r : A lto g e t h e r
S tric tly sp e a k in g , th e m a p p in g Ih d o e s n o t p ro d u c e a ( p , n ) -f-s e t w h ic h e x a c tly re p re s e n ts x h = y ( h , x0). In ste a d , it re tu rn s so m e b ig g e r se t [xh] su ch th a t x h is c o n ta in e d in it. O f c o u rse , w e a re in te re s te d in o b ta in in g a re s u lt as c lo s e as p o s s ib le to
F
0
C71
u<1 Springer
Found Com put M ath (2018) 18:1299-1332 1311
th e se t o f tru e so lu tio n s re p re s e n te d b y p ( h , x0). S o, fo r te c h n ic a l re a s o n s w h ic h w ill b e a p p a re n t in S ect. 2 .3 , w e d e c o m p o s e I h in to I h = @ + R su ch th a t @ : R m ^ R m a n d R : P (R m) ^ P (R m ). L e t @ ( x ) = $ a n d pu t:
$ i,[k] : = x i — 1,[k], i e {2,
$ 1-[0] := x 0,[0],
4> 1-[k] ; = I . F [k—1] (x p, [0]
k V
n
$ 0,[0]
s
w
IIk= 0
, x p ,[k—1] $ 1,[0] • 1,[k —1]
, k e {1, (9) (1 0) , n},
(1 1) (1 2)
L e t R ( x ) = r a n d pu t:
r i,[n+1] r 1, [p+1] r0,[0]
x i —1, [n+1]
= a
i e {2, [0, h ],
, p } ,
= h ■ (a* + b* ■ [0, h ] ) .
T h e m a p @ is c a lle d th e T aylo r p a r t , w h ile th e m a p R is c a lle d the R e m a i n d e r p a r t . T h is d e c o m p o s itio n is im p o rta n t fo r an effic ie n t re d u c tio n o f n e g a tiv e e ffe c ts c a u se d b y u s in g in te rv a l a rith m e tic , p rim a rily th e w ra p p in g effect, b u t also th e d e p e n d e n c y p ro b le m to so m e e x ten t.
b*
2 .3 R e d u c i n g t h e W r a p p i n g E f f e c t
A re p re s e n ta tio n o f o b je c ts in c o m p u ta tio n s as th e in te rv a l sets h a s its d ra w b a c k s.
P o s s ib ly th e m o s t se v e re o f th e m are th e p h e n o m e n a c a lle d th e w r a p p i n g effect a n d the d e p e n d e n c y p r o b l e m . T h e ir in flu e n c e is so d o m in a n t th a t th e y a re d is c u s s e d in v irtu a lly e v ery p a p e r in th e field o f rig o ro u s c o m p u ta tio n s (see [3 5 ] a n d re fe re n c e s th e re in ).
T h e d e p e n d e n c y p ro b le m a rises in in te rv a l a rith m e tic w h e n tw o v a lu e s th e o re tic a lly re p re s e n tin g th e sa m e (o r d e p e n d e n t) v a lu e are c o m b in e d . T h e m o s t triv ia l e x a m p le is an o p e ra tio n x x w h ic h is a lw a y s 0, b u t it is n o t th e c a s e fo r in te rv a ls. F o r e x a m p le , a p p ly in g th e o p e ra tio n to th e in te rv a l x = [ — 1, 1] g iv es as th e re s u lt the in te rv a l [ — 2, 2] w h ic h c o n ta in s 0, b u t it is far b ig g e r th a n w e w o u ld lik e it to be.
T h e w ra p p in g e ffe c t arise s w h e n o n e in te n d s to re p r e s e n t a re s u lt o f so m e e v alu a tio n on sets as a sim p le in te rv a l set. F ig u re 2 illu s tra te s th is w h e n w e c o n s id e r th e ro ta tio n o f th e sq u are.
O n e o f th e m o s tly u s e d an d e ffic ie n t m e th o d s fo r re d u c in g th e im p a c t o f th e w ra p p in g e ffe c t a n d th e d e p e n d e n c y p ro b le m w as p ro p o s e d b y L o h n e r [ 15]. In th e c o n te x t o f th e ite ra tio n o f m a p s a n d th e in te g ra tio n o f O D E s, h e p ro p o s e d to re p re s e n t sets b y p a ra lle lo g ra m s , i.e., in te rv a l sets in o th e r c o o rd in a te sy ste m s. In th e seq u el, w e fo llo w [3 5 ], an d w e sk e tc h th e L o h n e r m e th o d s briefly.
B y J w e d e n o te a c o m p u ta tio n o f I h (■) u s in g p o in t-w is e e v a lu a tio n o f th e T aylor p a rt, i.e.,
FoCTI u
< 1 S p rin g e r
Fig. 2 An illustration of the wrapping effect problem for a classical, idealized mathematical pendulum ODE x = —x . The picture shows a set of solutions in a phase space (x, x ). The gray boxes present points of initial box moved by the flow. The colored boxes present the wrapping effect occurring at each step when we want to enclose the moving points in a product of intervals in the basic coordinate system. For example, the blue square on the left encloses the image of the first iteration. Its image is then presented with blue rhombus which is enclosed again by an orange square. Then the process goes on (Color figure online)
J ([ x ]) : = I U ® ( x ) I + R ( [ x ]).
\xe[x] )
L e t u s c o n s id e r an ite ra tio n :
[xk] = J ([xk—1] ) , k e N
w ith in itia l set [x0].
L e t u s d e n o te x k = m i d ( [ x k ]) a n d [rk] = [xk] — x k . B y a sim p le a rg u m e n t b a se d on th e m e a n v a lu e th e o re m [3 5 ], it c a n b e sh o w n that:
[x k +1] C @ ( x k ) + [D 0 ([xk])] ■ [rk] + R ( [ x k ]).
W e c a n re f o rm u la te th e p ro b le m o f c o m p u tin g [xk+ i ] to th e fo llo w in g sy s te m o f e q u a tio n s:
[A k] = [ D 0 ( [ X k ] ) ] , (13)
xk+ 1 = m i d ( [ 0 ( x k ) + R ( [ X k ] ) ] ) , (14)
[Zk+1] = @([xk ]) + R ( [ x k ]) — xk+1, (15)
[rk+1] = [A k] ■ [rk] + [Zk+1]. (16)
N o w th e re d u c tio n in th e w ra p p in g e ffe c t c o u ld b e o b ta in e d b y c h o o sin g su ita b le re p re s e n ta tio n s o f sets [rk] an d a c a re fu l e v a lu a tio n o f E q . ( 16) . T h e te rm in o lo g y u se d fo r th is in [3 5 ] is t he r e a r r a n g e m e n t c o mp u t a t i o n s . W e w ill b rie fly d is c u ss p o s s ib le m e th o d s o f h a n d lin g E q . ( 16) .
F0C71
<1 Springer
u
Found Com put M ath (2018) 18:1299-1332 1313
M e t h o d 0 (In te rv a l S et): R e p re s e n ta tio n o f [rk] b y an in te rv a l b o x a n d th e d ire c t e v a lu a tio n o f ( 16) is e q u iv a le n t to d ire c tly c o m p u tin g I h ([xk+ 1]). T h is m e th o d is c a lle d a n in te r v a l s e t an d is th e le a s t effectiv e.
M e t h o d 1 (P a ra lle le p ip e d ): w e re q u ire th a t [rk] = B k ■ [rk ] fo r [rk ] b e in g an in terv al b o x a n d B k b e in g an in v e rtib le m a trix . T h e n ( 16) b e c o m e s:
[rk+1] = Bk+1 ^ B—+ 1 [A k ] Bk+1 ' [rk ] + B—_ \[zk+1]) .
S in c e it is d iffic u lt to o b ta in th e e x a c t m a trix in v e rse in c o m p u te r c a lc u la tio n s, w e w ill u s e in te rv a l m a tric e s [B k ] a n d [B—1] th a t c o n ta in B k an d B—1, re sp e c tiv e ly . T h u s th e e q u a tio n on [r] b e c o m e s:
[r:k+1] = ( [ Bk+ \] [ A k ] [ B k + 1 ^ ■ [rk] + [B k +1] 1[Zk+1].
I f [B k ] ’s a re w e ll-c h o se n , th e n th e fo rm u la in b ra c k e ts c a n b e e v a lu a te d to p ro d u c e a m a trix c lo s e to id e n tity w ith v e ry sm a ll d ia m e te r, th u s th e w ra p p in g e ffe c t re d u c tio n is a c h ie v e d . T h e P a ra lle le p ip e d m e th o d is o b ta in e d w h e n B k+ 1 is c h o se n su ch th at Bk+ 1 e [A k] [ Bk+ 1]. T h is a p p ro a c h is o f lim ite d u s e b e c a u s e o f th e n e e d to c o m p u te th e m a trix in v e rse o f a g e n e ra l m a trix B k+1, w h ic h m a y fail o r p ro d u c e u n a c c e p ta b le re s u lts i f Bk+ 1 is ill-c o n d itio n e d .
M e t h o d 2 (C u b o id ): th is is a m o d ific a tio n o f M e th o d 1. In th is m e th o d , w e c h o o se U e [A k ] [ Bk+ 1], an d w e d o th e flo a tin g p o in t a p p ro x im a te Q R d e c o m p o s itio n o f U = Q ■ R , w h e re Q is c lo s e to an o rth o g o n a l m a trix . N e x t w e o b ta in m a trix [ Q ] b y a p p ly in g th e in te rv a l (rig o ro u s ) G r a m - S c h m id t m e th o d to Q , so th e re ex is t o rth o g o n a l m a trix Q e [Q ] a n d Q —1 = Q T e [Q ] T .W e se t B k +1 = [Q ], B— = [ Q ] T .
M e t h o d 3 (D o u b le to n ): th is re p re s e n ta tio n is u s e d in o u r c o m p u ta tio n s as it p ro v e d to b e th e m o s t effic ie n t in n u m e ric a l tests [2 4 ] a n d in o th e r a p p lic a tio n s; se e [3 5 ] and re fe re n c e s th e re in . T h e o rig in a l id e a b y L o h n e r is to s e p a ra te th e e rro rs in tro d u c e d d u e to th e la rg e siz e o f in itia l d a ta an d th e lo c a l e rro rs in tro d u c e d b y th e n u m e ric a l m e th o d a t e v e ry step. N a m e ly , w e set:
[rk +1] = [E k + 1 ][r0 ] + [ k +1] [rk +1] = [Ak ][rk ] + [Zk+1] [E k +1] = [Ak ][ E k ] E0 = I d
w h e re [rk ] is e v a lu a te d b y an y m e th o d 0 - 2 . To re d u c e th e p o s s ib le w ra p p in g e ffe c t in th e p ro d u c t [A k ] [ E k ], L o h n e r p ro p o s e d th e fo llo w in g :
[rk +1] = C k+1[r0] + [ k +1]
[ ^ +1] = [A k ][rk ] + [Zk+1] + ([A k]C k — C k +1) [r0] r0 = 0 C0 = I d C k +1 e [A k ]C k .
A g a in , [I1] is e v a lu a te d b y an y m e th o d 0 - 2 . P le a s e n o te th a t th e re is n o n e e d to in v e rse a m a trix in th e d o u b le to n re p re s e n ta tio n w h e n [I1] is e v a lu a te d e ith e r b y M e th o d 0 or
FoCTI u
<1 Springer
M e th o d 2, so th is a p p ro a c h is su ita b le fo r th e c a s e w h e re A k m a y b e c lo s e to sin g u lar.
In th e c o m p u te r-a ssiste d p ro o fs p re s e n te d in th is p a p e r, w e u s e M e th o d 0 to re p re s e n t [r] b e c a u s e it is less c o m p u ta tio n a lly e x p e n siv e an d , in o u r c u rre n t se ttin g , u s in g th e o th e r m e th o d s h a v e n o t im p ro v e d th e re s u lts . T h is is p u z z lin g , as it c o n tra d ic ts our e x p e rie n c e w ith O D E s w h e re M e th o d 2 is p re fe ra b le , an d th u s it m ig h t b e w o rth w h ile to in v e s tig a te th is p h e n o m e n a in so m e la te r study.
2 .4 O p tim iz a t io n E x p lo iti n g t h e B lo c k S t r u c t u r e o f D ® ( x )
In o u r se ttin g , [A k ] = D @ ( x k ), w h e re x k is th e (p , n ) - f - s e t in th e k t h step o f in te g ra tio n . A s [A k ] is m x m m a trix , w h e re m is th e siz e o f (p , n )- re p re se n ta tio n (m = p ■ (n + 2) + 1), a n d th e re fo re , i f w e d e c id e to re p re s e n t th e e rro r p a rt [r] in d o u b le to n b y in te rv a l b o x (L o h n e r M e th o d 0 ), th e n th e m a trix m u ltip lic a tio n s in v o lv in g th e m a trix [A k ] ta k e th e m o s t o f th e e x e c u tio n tim e in o n e step o f in te g ra tio n in th e L o h n e r a lg o rith m , e s p e c ia lly fo r la rg e n a n d /o r p . F ro m E q s. (9- 12) , w e see th a t D @ ( x ) h a s a n ic e b lo c k stru c tu re an d c o n ta in s a la rg e n u m b e r o f ze ro e n trie s, i.e., it is a sp a rs e m a trix . T h is stru c tu re is w e ll v is ib le w h e n w e u s e th e fo llo w in g in d e x fu n c tio n I:
I( 0 , 0) = 1,
I ( i , k ) = 1 + (i — 1) ■ (n + 1) + k , 1 < i < p , 0 < k < n, I ( i , n + 1) = 1 + p ■ (n + 1) + i, 1 < i < p .
U n d e r th is in d e x fu n c tio n , [A k ] is o f th e form :
A11o 0 A 13 0
iiS
I d 0 0
0 0 0 0
w h e re A11 is n + 1 x 1 m a trix (c o lu m n v e c to r), A13 is n + 1 x n + 2, a n d I d is an id e n tity m a trix o f siz e (p — 1) ■ (n + 1 ) . T h e re fo re , w e u s e th is in d e x fu n c tio n to d efin e b lo c k s fo r all m a tric e s an d v e c to rs a p p e a rin g in all m e th o d s d is c u s s e d in S ect. 2 .3 , as, w ith su c h b lo c k re p re s e n ta tio n o f m a tric e s an d v e c to rs, w e c a n e a s ily p ro g ra m th e m u ltip lic a tio n b y a m a trix or a v e c to r so th a t all th e o p e ra tio n s on an y z e ro b lo c k are a v o id e d . W e w ill re fe r to th is as th e o p tim iz e d a lg o rith m .
I f w e h a v e an a rb itra ry m a trix C , th e n th e c o s t o f c o m p u tin g [A k] ■ C b y a s ta n d a rd a lg o rith m fo r th e m a trix m u ltip lic a tio n is o f o rd e r O ( n 3 ■ p 3) in b o th th e scalar a d d itio n an d m u ltip lic a tio n o p e ra tio n s (w e re m in d , th a t p , n a re th e p a ra m e te rs o f ( p , n )- re p re s e n ta tio n ). In th e c a s e o f th e o p tim iz e d a lg o rith m , th e b lo c k stru c tu re and sp a rse n e ss o f [A k] re d u c e th e c o m p u ta tio n a l c o s t to O ( n 2 ■ p 2) in sc a la r a d d itio n s and 0 ( n 3) in sc a la r m u ltip lic a tio n s.
T h e c o m p u ta tio n tim e s fo r th e c o m p u te r-a s s is te d p ro o fs d is c u s s e d in S ect. 4 on the 2 .5 0 G H z p ro c e s s o r (se e S ect. 4 fo r a d e ta ile d sp e c ific a tio n ) a re p re s e n te d in T ab le 1.
W e se e th a t th e o p tim iz e d a lg o rith m is m u c h fa s te r th an th e d ire c t m u ltip lic a tio n , th e sp e e d -u p is e v id e n t e s p e c ia lly fo r th e la rg e r (p , n )-re p re s e n ta tio n s .
FoCTI
1
Springeru
Found Com put M ath (2018) 18:1299-1332 1315
Table 1 A comparison of the execution times for computer-assisted proofs for standard and optimized matrix multiplication on 2.5 GHz processor (no multi-threading)
Proof (p , p) Standard multiplication Optimized multiplication
Theorem 12 (32, 4) 45 s 1 2s
Theorem 13 (128, 4) 133 min 12 min
3 Poincare Map for Delay Differential Equations
3 .1 D e fin itio n o f a P o in c a r e M a p
W e b e g in w ith th e d e fin itio n o f th e (tra n sv e rsa l) se c tio n o f th e sem iflo w p a ss o c ia te d w ith ( 1). F irst, w e w o u ld lik e to re c a ll th e O D E se ttin g w h e re , fo r a flow p : R x R m ^ R m , a (lo c a l) tra n sv e rs a l se c tio n S is u s u a lly d e fin e d as a (s u b se t o f) sm o o th m a n ifo ld M o f c o d im e n sio n o n e sa tisfy in g th e tra n s v e rs a lity co n d itio n :
d
— p ( t , x0 ) |t= 0 = f ( x0 ) e Tx0M , x0 e S, (17)
w h e re Tx M d e n o te s th e ta n g e n t b u n d le at x . I f M is a h y p e rp la n e
M = {x e R m | s ■ x = a}
fo r so m e g iv en n o rm a l v e c to r s e R m , ||s || = 0 a n d a e R , th en c o n d itio n ( 17) b e c o m e s
s ■ f ( x0 ) = 0, x0 e S.
W e w ill u s e a sim ila r a p p ro a c h in th e c o n te x t o f th e sem iflo w p a ss o c ia te d w ith ( 1) . W e w ill re s tric t o u rs e lv e s to th e lin e a r se c tio n s, an d w e w ill u s e th e fa c t th a t fro m E q . (4 ) an d L e m m a 2 , it fo llo w s
d
d t p ( t , x0) | t= t0 x t0
fo r an y t0 > t . M o re o v e r, xt0 is o f cla ss Cp—1 w h e re v e r t0 > n ■ t . T h is o b se rv a tio n w ill b e c ru c ia l fo r th e d e fin itio n o f a tra n sv e rs a l se c tio n in th e D D E c o n te x t and, later, fo r th e rig o ro u s c o m p u ta tio n o f P o in c a re m ap s.
D e f in itio n 6 L e t p b e th e sem iflo w a ss o c ia te d w ith th e sy s te m ( 1) . L e t n e N , s : Cp ^ R b e a c o n tin u o u s affin e m a p p in g , i.e ., s ( x ) = l ( x ) + a, w h e re l is a b o u n d e d lin e a r fu n c tio n a l a n d a e R . W e d e fin e a g l o b a l C n -s e c tio n as a h y p e rp la n e :
S = {x e Cp | s ( x ) = 0}.
A n y c o n v e x a n d b o u n d e d su b s e t S c S is c a lle d a l o c a l C n -s e c tio n (o r sim p ly a se c tio n ).
FoCTI u
<1 Springer
A se c tio n S is sa id to b e t r a n s v e r s a l i f th e re e x is ts a c o n v e x o p en set W c C n such th a t
W n 5 = U, W = W —U U U W +, w h e re
W — = (x e W | s ( x ) < 0}, W + = (x e W | s ( x ) > 0}, c l ( U ) = cl ( S ),
s a tisfy in g th e c o n d itio n
l ( x ) > 0, Vx e W n C n+ 1. (18)
W e w ill re f e r to ( 18) as th e tr a n s v e r s a lit y conditio n.
R e m a r k 4 P le a s e n o te th a t th e re q u ire m e n t x e W n Cn + 1 in ( 18) is e s s e n tia l to g u a ra n te e th a t x a n d th u s l ( x ) in ( 18) a re w ell d efin ed , as, fo r x0 e Cn + 1 a n d t > 0, it m ig h t h a p p e n th a t y ( t , x 0) is o f c la s s C k , k < n + 1 (th e lo ss o f re g u la rity ), but L e m m a 2 states th a t w e o n ly n e e d to “lo n g e n o u g h ” in te g ra te th e in itia l fu n c tio n s to get rid o f th is p ro b le m co m p le te ly . T h e s e tw o p h e n o m e n a a re illu s tra te d in th e fo llo w in g e x am p le.
L e t x0 e W fo r W as in D e fin itio n 6. In g e n e ra l, it h a p p e n s th a t y ( e , x 0) / C n for sm a ll e > 0. T h is m a y s e e m at first to c o n tra d ic t in tu itio n fro m L e m m a 2 , b u t in fact it is n o t. C o n s id e r th e fo llo w in g r.h .s. o f E q . ( 1) :
f (Z1 , Z2 ) = Z1, V z1, Z2 e R
L e t x0 = 1 b e an in itia l fu n c tio n a n d le t x b e a so lu tio n to E q . ( 1) w ith x |[—1,0] = x0
a n d d e la y t = 1. W e see th a t x ( t ) is C “ on [—1, 0 ). H o w ev er, at t = 0, w e h a v e
0 = x ' (0—) = 1 = x ' (0+ ) = f ( x ( — 1), x ( 0 ) ) ,
so x t : [ - 1 , 0] —— R is o n ly o f class C0 fo r an y t e (0, 1). T h is is a v e ry u n d e sira b le p h e n o m e n a , b u t th e so lu tio n w ill b e s m o o th e d a fte r a fu ll d elay , a c c o rd in g to L e m m a 2 . A s x ( t ) = 1 + 1 on [0, 1], w e h a v e at t = 1:
1 = x ' ( r ) = x ' (1+ ) = f ( x (0), x (1) ) .
O n e c a n sh o w ag a in th a t x (2)(1 —) = x (2) (1 + ), th e re fo re x is o f c la s s C1 on (1, 2 ), an d th e sm o o th in g o f so lu tio n s g o e s on w ith th e in c re a s in g t .
T h is sh o w s fo r a n y x0 e C n , i f m > n ■ t, th e n w e h a v e “ o n ly ” y ( m , x 0) e C n in a g e n e ra l ca se . O n th e o th e r h a n d , “ lo n g e n o u g h ” in te g ra tio n tim e m c a n b e u se d to g u a ra n te e th a t e v e ry in itia l fu n c tio n x e C0 h a s a w e ll d efin ed im a g e in C n u n d e r m a p p in g y ( m , ■). T h is is e s s e n tia l in th e fo llo w in g c o n s tru c tio n o f a P o in c a re m ap for D D E s (T h e o re m 5 a n d D e fin itio n 7 ) .
T h e o r e m 5 A s s u m e n e N, V c C 0. L e t S be a lo c a l tr a n s v e r s a l C n - s e c t i o n f o r ( 1) a n d let W b e a s in D e fin i ti o n 6 . L e t m = (n + 1) ■ t.
FoCJI
<1 Springer
u
Found Com put M ath (2018) 18:1299-1332 1317
A s s u m e th a t th ere e x is t t 1, t2 e R , a < t1 < t2, su c h th a t th e fo llo w in g c o n d itio n s h o ld f o r a ll x e V :
p ( [ t 1, t2], x ) C W , <p(t1, x ) e W — a n d p ( t 2 , x ) e W + . (19)
Then, f o r ea ch x 0 e V , th ere e x is ts u n iq u e t S ( x 0) e ( t1, t2) su c h th a t p (tS ( x 0), x 0) e S, A lso , tS : V ^ [ t 1, t2] is co n tin u o u s,
P r o o f L e t x0 e V . B y a ss u m p tio n s, p ( t , x 0) = x t e W fo r t e [t1, t2] b u t also x t e C n+ 1, b y th e a ss u m p tio n on c o n sta n ts a , t 1, t2 (b y L e m m a 2 ). So s ( x t ) is w ell d e fin e d an d C o n d itio n ( 18) g u a ra n te e s th a t
d ( d \
— s ( p ( t , x )) = l ( — p ( t , x ) \ = l ( x t ) > 0, t e [ t1, t2].
T h e re fo re , th e fu n c tio n d e fin e d b y s ( t ) : = s ( x t ) is c o n tin u o u s a n d stric tly in c re a s in g on [ t1, t2]. N ow , fro m ( 19), it fo llo w s th e re e x ists u n iq u e t0 e ( t 1, t2) su c h th at s ( p ( t 0 , x 0)) = c. T o g e th e r w ith c o n tin u ity o f p (L e m m a 1) , th is im p lie s c o n tin u ity o f
t s : V ^ ( t1, t2 ). □
D e f in itio n 7 T h e s a m e a ss u m p tio n s as in T h e o re m 5 , in p a rtic u la r a ss u m e t1 > t ■ (n + 1) = a . W e d e fin e th e tr a n s itio n m a p to th e s e c tio n S a fte r (a t lea st) a b y
P>a : V ^ S C C n , P > a (x0 ) : = p ( t s ( x 0 ), x 0 ) ,
fo r tS d e fin e d as in T h e o re m 5 . I f V c S, th e n th e m a p P>a w ill b e c a lle d th e P o in ca re re tu rn m a p o n th e s e c tio n S a fte r a .
F in a lly , w e s ta te th e la s t an d th e m o s t im p o rta n t th e o re m th a t w ill allo w u s to a p p ly to p o lo g ic a l fix ed p o in t th e o re m s to P>a .
T h e o r e m 6 C o n s id e r P o in ca re m a p (a fte r a ) P>a : S D V ^ S f o r so m e s e c tio n S u n d e r th e sa m e a ss u m p tio n s a s in T h e o re m 5, e s p e c ia lly a ss u m e (n + 1) ■ t < a <
t S ( V ) e [t1, t2],
A s s u m e a d d itio n a lly th a t p ( [ a , t2], V ) is b o u n d e d in C n ,
T h en th e m a p P>a is c o n tin u o u s a n d c o m p a c t in C n , i,e,, i f K c V is b o u n d ed , th e n c l ( P>a ( K )) is c o m p a c t in C n ,
P r o o f B y T h e o re m 5 , P>a is w e ll d e fin e d fo r a n y x0 e V sin c e t1 > a > (n + 1) ■ t
a n d P>a (x0) e C n+ 1.
T h e c o n tin u ity fo llo w s im m e d ia te ly fro m th e c o n tin u ity o f tS (T h e o re m 5 ) an d p (L e m m a 1) .
L e t D = P>a ( V ). F ro m o u r a ss u m p tio n s, it fo llo w s th a t D is b o u n d e d in C n . A k n o w n c o n s e q u e n c e o f th e A r z e la - A s c o li T h e o re m is th at, i f D C C n is c lo s e d an d b o u n d e d , x e D , x (n+ 1) e x ists, an d th e re is M su ch th a t su p t |x (n+ 1) ( t ) < M fo r all x e D , th e n D is c o m p a c t (in C n -n o rm ). T h e re fo re , to fin ish th e p ro o f, it is e n o u g h to sh o w th a t th e re is a u n ifo rm b o u n d on P>a ( x ) (n+ 1) . F o r th is, it is su ffic ie n t to h a v e a u n if o rm b o u n d on ( p ( t , x ) ) (n+ 1) fo r t e [ a , su p x e V t S ]. T h e e x is te n c e o f th is b o u n d fo llo w s fro m b o u n d e d n e s s o f d eriv ativ es u p to o rd e r n a n d fo rm u la (7 ) . □
FoCTI u
< 1 S p rin g e r
