Proper Orthogonal Decomposition and Dynamic Mode Decomposition of Delay-Differential Equations

Proceedings, 14th on Delay Systems Pesti Vigadó, Budapest, Hungary, June 28-30, 2018 Proceedings, 14th IFAC IFAC Workshop Workshop on Time Time Delay Systems Proceedings, 14th IFAC Workshop on Time Delay Systems Pesti Vigadó, Vigadó, Budapest, Hungary, 28-30, 2018 Pesti Budapest, Hungary, June June 28-30, 2018 Available online at www.sciencedirect.com Pesti Vigadó, Budapest, Hungary, June 28-30, 2018 Proceedings, 14th IFAC Workshop on Time Delay Systems Pesti Vigadó, Budapest, Hungary, June 28-30, 2018

ScienceDirect

IFAC PapersOnLine 51-14 (2018) 254–258 Proper Orthogonal Decomposition and Proper Orthogonal Decomposition and Proper Orthogonal Decompositionofand Dynamic Mode Decomposition Dynamic Mode Decomposition Proper Orthogonal Decompositionof and Dynamic Mode Decomposition Delay-Differential Equations of Delay-Differential Equations of Dynamic Mode Decomposition Delay-Differential Equations ∗ B. Heizer, T. Kalmar-Nagy Delay-Differential Equations B. Heizer, T. Kalmar-Nagy ∗∗

B. Heizer, T. Kalmar-Nagy ∗ Kalmar-Nagy Faculty of Mechanical Engineering, ∗ Kalmar-Nagy Faculty of Engineering, Faculty of Mechanical Mechanical Engineering, and Economics, Budapest, Hungary Department of Fluid Mechanics, Faculty of Mechanical Engineering, Budapest University of Technology and Economics, Budapest, Budapest University of Technology and Economics, Budapest, Hungary Hungary (e-mail: [email protected]) ∗ Budapest University of Mechanics, Technology and Economics, Budapest, Hungary Department of Fluid Faculty of Mechanical Engineering, (e-mail: [email protected]) (e-mail: [email protected]) (e-mail: [email protected]) Budapest University of Technology and Economics, Budapest, Hungary [email protected]) Abstract: Our work is based (e-mail: on the observation that delay-differential equations can be recast Abstract: Our work work is isequations based on on the the observation that delay-differential be Abstract: Our based observation that delay-differential equations can be recast recast as partial differential of two variables. Thus the numericalequations solution can of the delay Abstract: Our work isequations based on the observation that delay-differential equations can be recast as partialcan differential ofspatio-temporal two variables. Thus theand numerical solution of the the delay as partial differential equations two variables. Thus the numerical of delay equation be thought of as aof process, snapshotsolution based methods like as partialcan differential equations ofspatio-temporal two variables. Thus theand numerical solution of the delay Abstract: Our work is based on the observation that delay-differential equations can be recast equation be thought of as a process, snapshot based methods like equation can be thought of as a spatio-temporal process, and snapshot Proper Orthogonal Decomposition and Dynamic Mode Decomposition allowbased us to methods compute like the equation can be thought of as aofspatio-temporal process, snapshot based methods like as partial differential equations twoDynamic variables. Thus theand numerical solution of the delay Proper Orthogonal Decomposition and Dynamic Mode Decomposition allow us to compute the Proper Orthogonal Decomposition and Mode Decomposition allow us to compute eigenvalues of the evolution operator. The rightmost eigenvalues of the Hayes equation canthe be Proper Orthogonal Decomposition andThe Dynamic Mode Decomposition usequation to methods compute the equation can of operator. as a spatio-temporal process, and snapshot based eigenvalues of be thethought evolution operator. rightmost eigenvalues of the theallow Hayes canlike be eigenvalues of the evolution The rightmost eigenvalues of Hayes equation can be well approximated. eigenvalues of the evolution operator. The rightmost eigenvalues of theallow Hayes equation canthe be Proper Orthogonal Decomposition and Dynamic Mode Decomposition us to compute well approximated. well approximated. well approximated. eigenvalues the evolution operator. rightmost eigenvalues the Hayes equation can be © 2018, IFACof(International Federation of The Automatic Control) Hosting byofElsevier Ltd. All rights reserved. Keywords: Delay-differential equations, Proper Orthogonal Decomposition, Dynamic Mode well approximated. Keywords: equations, Keywords: Delay-differential Delay-differential equations, Proper Proper Orthogonal Orthogonal Decomposition, Decomposition, Dynamic Dynamic Mode Mode Decomposition, Rightmost eigenvalues. Keywords: Delay-differential equations, Proper Orthogonal Decomposition, Dynamic Mode Decomposition, Rightmost eigenvalues. Decomposition, Rightmost eigenvalues. Decomposition, Rightmost eigenvalues. Keywords: Delay-differential equations, Proper Orthogonal Decomposition, Dynamic Mode 2. PROPER ORTHOGONAL DECOMPOSITION 1. INTRODUCTION Decomposition, Rightmost eigenvalues. 2. 1. INTRODUCTION INTRODUCTION 2. PROPER PROPER ORTHOGONAL ORTHOGONAL DECOMPOSITION DECOMPOSITION 1. 2. PROPER ORTHOGONAL DECOMPOSITION 1. INTRODUCTION Analysis of delay-differential equations is an ever-growing Proper orthogonal decomposition (POD) is a powerful 2. PROPER ORTHOGONAL DECOMPOSITION 1. INTRODUCTION Analysis ofgoal delay-differential equations is an an ever-growing decomposition is Analysis delay-differential equations is ever-growing field. Ourof is to understand whether delay equations Proper Proper orthogonal orthogonal decomposition (POD) is a a powerful powerful dimensionality reduction technique.(POD) The singular value Analysis of delay-differential equations is an ever-growing orthogonal decomposition (POD) is a powerful field. Our goal understood is to to understand understand whether delay delay equations dimensionality reduction technique. The singular value field. Our goal is whether equations can be better or represented in terms of a Proper dimensionality reduction technique. The singular value decomposition field. Our is to understand whether delay equations Analysis ofgoal delay-differential is an dimensionality reduction technique. The singular value orthogonal decomposition (POD) is a powerful can be better understood orequations represented inever-growing terms of aa Proper decomposition can be better understood or represented in of natural “basis” (Asl and Ulsoy [2003], Amann etterms al. [2007], ∗ decomposition X = U ΣV (8) can be better understood or [2003], represented in etterms of a dimensionality field. Our goal is to understand whether delay equations decomposition reduction technique. The singular value natural “basis” (Asl and Ulsoy Amann al. [2007], ∗ natural (Asl andThe Ulsoy [2003],form Amann et al. [2007], ∗ Michiels“basis” et al. [2011]). general of scalar delayX = U ΣV (8) X = U ΣV (8) natural “basis” (Asl and Ulsoy [2003], Amann et al. [2007], ∗ can be better understood or represented in terms of a Michiels et equations al. [2011]). [2011]).is The The general general form form of of scalar scalar delaydelay- decomposition yields the matrix U that X =contains U ΣV the orthogonal spatial (8) Michiels et al. differential yields the matrix U that contains the orthogonal spatial Michiels et al. [2011]). The general form of scalar delaynatural “basis” (Asl and Ulsoy [2003], Amann et al. [2007], ∗ yields matrix U that orthogonal differential equations equations is is modes the (POD basis), ΣXthat the singular spatial values =contains Ucontains ΣV the (8) differential dx(t) yields the matrix U that contains the orthogonal spatial modes (POD basis), Σ that that information contains the the issingular singular values differential is=The Michiels et equations al. [2011]). general form of scalar delaymodes (POD basis), Σ contains values σ , while all the temporal contained in f (x(t), x(t − τ )), (1) i dx(t) dx(t) modes (POD basis), Σ that contains singular values the matrix that contains the the orthogonal spatial ∗ U dt is= σ all the temporal information is contained in = ff (x(t), (x(t), x(t x(t − − ττ )), )), (1) yields differential equations ii ,, while dx(t) σ while all the temporal information is contained in matrix V and ( denotes the conjugate transpose). The (1) ∗ dt = = σi , while all the temporal information issingular contained in modes (POD basis), Σ that contains the values f (x(t), x(t − τ )), (1) ∗ denotes dt matrix V and ( the conjugate transpose). The x(t) θ(t), −τ ≤ t ≤ 0, (2) matrix Vvalues and (in denotes theintuitive conjugate transpose). The dx(t) singular Σ have an meaning, they are ∗ dt x(t) == θ(t), −τ ≤ tt ≤ 0, (2) matrix Vvalues andthe (in denotes theinformation conjugate transpose). The σ , while all temporal is contained in f (x(t), − τintroducing )), (1) i x(t) = θ(t), −τx(t ≤By ≤ 0, (2) singular Σ have an intuitive meaning, they are where θ(t) is the initial function. the sosingular values in have an intuitive are of the proportional with the energy content meaning, (L1 normthey ∗ Σ x(t) = θ(t), −τ ≤By t ≤introducing 0, (2) dt singular values in Σ have an intuitive meaning, they are matrix V and ( denotes the conjugate transpose). The where θ(t) is the initial function. the sothe proportional with energy (L where θ(t) initial function. the so- singular called shift isofthe time u(t, s) =−τ x(t≤By +t s), value norm ofWith the proportional with the energy content contentmodes (L11 norm values) of the the corresponding in U .of x(t) = θ(t), ≤introducing 0,the initial (2) where θ(t) isofthe initial function. By introducing the so- singular ofWith the proportional with energy contentmodes (L1 norm values inof Σthe have an intuitive meaning, they are called u(t, = + initial value values) the corresponding in U ..obtain called shift shift time u(t,bes) s)recast = x(t x(tinto + s), s), the initialpartial value singular problem Eqs.of(1,time 2) can the the following singular values) of the corresponding modes in U With the orthogonal basis contained in U , it is possible to called θ(t) shift time u(t,be s)recast = x(tinto + s), the initialthe value where isof(1, the initial function. By introducing so- singular values) of the corresponding modes in U . With norm of the proportional with the energy content (L problem Eqs. 2) can the following partial 1 the orthogonal basis contained in U , it is possible to obtain problem Eqs. (1, 2) can be recast into the following partial differential equation (Bellen and Maset [2000]) the orthogonal basis contained in U , it isthrough possible to obtain a reduced order model of the DDE Galerkin problem Eqs. 2) can into the[2000]) following called shift of(1,time u(t,bes)recast = x(tMaset + s), the initialpartial value singular the orthogonal basis contained in U , it ismodes possible values) ofmodel the corresponding in to U .obtain With differential equation (Bellen and aprojection. reduced order of the DDE through Galerkin differential equation (Bellen and Maset [2000]) a reduced order model of the DDE through Galerkin The solution u(t, s) can be expanded in terms ∂u(t, s) ∂u(t, s) differential equation (Bellen and Maset [2000]) problem Eqs. (1, 2) can be recast into the following partial a reduced order model of the DDE through Galerkin the orthogonal basis contained in U , it is possible to obtain The u(t, s ∈ [−τ, 0], (3) projection. s) ∂u(t, s) = ∂u(t, projection. The solution solution u(t, s) s) can can be be expanded expanded in in terms terms n POD modes as ∂u(t, s) ,, Maset ∂u(t, ∂s and ∂t s) differential equation [2000]) ∈ [−τ, [−τ, 0], (3) of = The s) can be expanded in terms aprojection. reduced ordersolution model u(t, ofn the DDE through Galerkin ∂u(t, s) (Bellen s) , ss ∈ = ∂u(t, 0], (3) of n POD modes as of n POD modes as ∂s ∂t  (3) of ∂t s) ∂s s) , s ∈ [−τ, 0], ∂u(t, n POD modes as projection. The solution u(t,nn s) can be expanded in terms ∂u(t, s) = ∂u(t, ∂t s) ai (t)φi (s), (9) u(t, s) ≈  =∂s f (u(t, 0), u(t, −τ0], )), (4)  ∂u(t,  , s ∈ [−τ, (3) = n ∂u(t, s)  of n POD modes as  ∂t a (9) u(t, s) ≈ f (u(t, 0), u(t, −τ )), (4) ii (t)φ ii (s), ∂s ∂t s) s=0 = ∂u(t, a (t)φ (s), (9) u(t, s) ≈ = f (u(t, 0), u(t, −τ )), (4) i=1 ∂t s=0 = f (u(t, 0), u(t, −τ )), n ai (t)φi (s), (9) u(t, s) ≈  (4) ∂t i=1 u(0, s) = θ(s), s ∈ [−τ, 0]. (5) s=0 ∂u(t, s) i=1 ∂t where n is the number of POD modes and the φ ’s belong  s=0 i a (t)φ (s), (9) u(t, s) ≈ = f (u(t, 0), u(t, −τ )), (4) u(0, s) = θ(s), s ∈ [−τ, 0]. (5) i=1 i i u(0, s) = θ(s), s ∈ [−τ, (5) where  can is POD modes φ The u(t, s) function be viewed as a0]. ∂t where northogonal is the the number number of POD POD modes and and the φii ’s ’s belong belong the n set ofof Thethe corresponding u(0, s)s=0 = θ(s), s ∈ [−τ, 0].spatio-temporal (5) to i=1 modes. The u(t, s) function can be viewed as a spatio-temporal where n is the number of POD modes and the φ ’s belong i The u(t, s) function be viewed as a0].spatio-temporal the orthogonal set modes. The corresponding representation of the value problem Eqs. (1, (5) 2), to toi (t)’s the can orthogonal set of of POD POD modes. The method. corresponding be determined by the Galerkin First, u(0, s)can =initial θ(s), s ∈ [−τ, The u(t, s) function can be viewed as a spatio-temporal representation of the the initial value problem Eqs.solution (1, 2), 2), a to the ncan orthogonal set ofof POD modes. Thethe corresponding where is the number POD modes and φi ’s belong representation of initial value problem Eqs. (1, a be determined by the Galerkin method. First, where the variable s is the spatial dimension. The ii (t)’s a (t)’s can be determined by the Galerkin method. First, substituting Eq. (9) into Eq. (3) we get representation of the initial value problem Eqs. (1, 2), The u(t, s) function can be viewed as a spatio-temporal where the variable is the the spatial dimension. The solution solution a beEq. determined by (3) the Galerkin First, the can orthogonal setinto of POD modes. The method. corresponding i (t)’s where ss is dimension. The substituting (9) Eq. we get of Eqs.the (1,variable 2) at time tk isspatial written into a column vector to substituting Eq. (9) into Eq. (3) we get n n where the variable s is the spatial dimension. The representation of time the initial value problem Eqs.solution (1, 2), a   of Eqs. 2) ttkk is into vector substituting Eq. (9) into Eq. we get∂φi (s) determined by (3) the Galerkin method. First, (t) ∂a i (t)’s can be of Eqs. a(1, (1, 2) at at time is written written into a a column column vector (called snapshot) i n n n φi (s) ∂a (t) =  n ai (t) ∂φ (s) . of Eqs. a(1, 2) at time tk isspatial written into a column vector substituting (10) where variable s is the dimension. The solution   (called i (t) i (s) Eq. (9) into Eq.=(3) n n we get ∂a∂t ∂φ∂s (calledthe a snapshot) snapshot) T   (10) φ a i (s) i (t) ∂aii (t) = ∂φii (s) .. x = at [u(ttime , s2 ), ...,into u(tk ,asncolumn )]T . vector (6) (10) φ (s) a (t) (called snapshot) of Eqs. a(1, tku(t is kwritten i=1 i=1 k2) k , s1 ), i i ∂t =  n φi (s) ∂t n ai (t) ∂s T . . (10) ∂s  x = [u(t , s ), u(t , s ), ..., u(t , s )] (6) i=1 i=1 k k 1 k 2 k n (t) (s) ∂a ∂φ [u(tk , sX s2 ),constructed ..., u(tk , sn )]Tfrom . (6) Projection of (called ax snapshot) i=1 i=1 i i k = 1 ), u(t k , be ∂t ∂s A snapshot matrix can these Eq. (10) onto=the φj , xk = [u(tk , s1 ), u(tk , s2 ), ..., u(tk , sn )] . (6) . (10) φi (s) ai (t) basis functions i=1 i=1POD A matrix X can constructed these Tfrom Projection Eq. onto POD ∂s functions A snapshot snapshot matrix X cank , be be constructed from these snapshots Projection of Eq. (10) (10) ∂t onto the the POD basis basis functions φ φj ,, j = 1, ..., n of yields x = [u(t , s ), u(t s ), ..., u(t , s )] . (6)   i=1 i=1 k k 1 2 k n A snapshot matrix X can be constructed from these Projection of Eq. (10) onto the POD basis functions φjj , snapshots j = 1, ..., n yields snapshots     j = 1, ..., n yields n  can be constructed   snapshots A snapshot matrix X from these Projection j= ..., n of yields Eq. 1(10) onto the POD φj , da1, (s) functions  ∂φbasis  j (t) n  x1 x2 . . . xm  . n ai (t)  ∂φi (s) , φj (s) . (11) X= (7)  = 1 da (t)   ∂φ∂s  snapshots j ii (s) n 1 da (t)  j = 1, ..., n yields  j x x . . . x   X = . (7)   φ dt (s), φ (s)  1 2 m = , φ a (t) (s) . (11) X =  x1 x2 . . . xm  . (7) daj (t) = j 1 j ai (t) ∂φi (s) , φjj (s) . (11) i=1 X = x1 x2 . . . xm . (7) dt = φ n aii (t)  ∂s , φj (s) . (11) dt φjj (s), (s),1φ φjj (s) (s)  ∂s i=1 da ∂φ∂s (t) i=1 j i (s) can be written φ dt (s), φ (s) This representation of the problem allows us to apply snap  j j differential equations X = x1 x2 . . . xm . (7) This set of = ordinary , φj (s) . (11) i=1 ai (t) This representation of the problem allows us to apply snapordinary equations This representation of the us to apply snap- This shot based approaches likeproblem Properallows Orthogonal Decompoφj (s), φj differential dtset ∂s can (s) i=1 This set of ofform ordinary differential equations can be be written written in matrix This representation of the problem allows us to apply snapThis set ofform ordinary differential equations can be written shot like Proper Decompomatrix shot based based approaches likeDynamic Proper Orthogonal Orthogonal Decompo- in sition (Kutzapproaches [2013]) and Mode Decomposition da in matrix form shot based approaches likeproblem Properallows Orthogonal DecompoThis ofand the usDecomposition to apply snap- This in matrix set ofform ordinary differential equations can be written sition (Kutz [2013]) Mode = La, (12) da sitionrepresentation (Kutz [2013]) Dynamic Mode (Schmid [2010], Kutz and et al.Dynamic [2016], Tu et al.Decomposition [2013]). da dt = = La, La, (12) sition (Kutz [2013]) and Dynamic Mode Decomposition shot based approaches like Proper Orthogonal Decompoda in matrix form (12) (Schmid [2010], Kutz et al. [2016], Tu et al. [2013]). (Schmid [2010], Kutz et al. [2016], Tu et al. [2013]). dt (12) dt = La, (Schmid [2010], Kutz and et al.Dynamic [2016], Tu et al.Decomposition [2013]). sition (Kutz [2013]) Mode da dt = La, (12) Copyright © 2018 IFAC 254 (Schmid [2010], Kutz et al. [2016], Tu et al. [2013]). 2405-8963 © 2018, IFAC (International Federation of Automatic Control) Hosting by Elsevier Ltd. All rights dt reserved. Copyright 254 Copyright © © 2018 2018 IFAC IFAC 254 B. Heizer, T. ∗ Department of Fluid Mechanics, ∗ B. Heizer, T. ∗ Department of Fluid of Fluid Mechanics, Budapest University of Mechanics, Technology ∗ Department

Peer review©under of International Federation of Automatic Copyright 2018 responsibility IFAC 254Control. 10.1016/j.ifacol.2018.07.232 Copyright © 2018 IFAC 254

2018 IFAC TDS Budapest, Hungary, June 28-30, 2018

Balázs Heizer et al. / IFAC PapersOnLine 51-14 (2018) 254–258

255

T

where a = [a1 , a2 , ..., an ] and the elements of L are given by   ∂φi (s) ∂s , φj (s) . (13) Lij = φj (s), φj (s) The boundary conditions (4,5) enter into matrix L through relationship (8) between POD modes φi (s) and the snapshot matrix X. The matrix L is a finite dimensional approximation of the infinitesimal generator of the DDE. By computing the spectra of L we can approximate the rightmost eigenvalues of the original system. 3. DYNAMIC MODE DECOMPOSITION Dynamic mode decomposition (DMD) provides an alternative way to construct reduced order model of dynamic system from numerical data. Initially DMD was developed to decompose fluid flows and find spatio-temporal coherent structures (Schmid [2010]), however it became popular soon in a more general context (Brunton et al. [2016], Hua et al. [2016], Proctor and Eckhoff [2015]). DMD is based on the assumption that the snapshots xk in the data matrix Eq. (7) are connected by a linear map (14) xk+1 = Axk . In order to compute matrix A we construct the two matrices   X 1 =  x1 x2 . . . xm−1  ,

(15)

X 2 =  x2 x3 . . . xm  .

(16)





from the data matrix X (see Eq. (7)). These matrices are connected by (assuming an underlying linear process) (17) X 2 = AX 1 . The singular value decomposition of X 1 is given by (18) X 1 = U 1 Σ1 V ∗1 . Matrix A is computed as (+ denotes the Moore-Penrose pseudoinverse) −1 ∗ (19) A = X 2X + 1 = X 2 V 1 Σ1 U 1 . If the underlying delay equation Eq. (1) is linear (i.e. f (x(t), x(t − τ )) = ax(t) + bx(t − τ )), then A is a finite dimensional approximation of its evolution operator. The eigenvalues of A approximate the spectrum of the evolution operator, which is related to the spectrum of Eq (1) by the exponential transform.

4. NUMERICAL RESULTS

Fig. 1. The normalized singular values of X φ(t) = 1. This numerical solution xnum (t) was used to generate POD and DMD results. 4.1 POD results The energy content of the POD modes are depicted in Fig. (1). The first four POD modes contain 98.87% of the energy, these are illustrated in Fig. (2). We used these four modes for the expansion (9) and after projecting Eq. (20) onto these modes we get the following system of ordinary differential equations (the initial function is similarly projected, resulting in the initial conditions)      −0.05 1.84 0.01 −0.01 a1 a1 d  a2   −3.24 −0.11 −0.02 −0.02   a2    =  −0.59 −3.24 −1.09 6.06   a  . (22) 3 dt a3 a4 a4 −1.75 −1.06 −10.51 −1.17 The solutions of Eq. (22) are shown Fig. (3). Eigenvalues of the coefficient matrix in Eq. (22) are approximating the eigenvalues of the evolution operator of Eq. (20). In Fig. (5) we compared these 4 approximate eigenvalues (2 complex pairs) with those obtained from the Matlab package published by Breda et al. [2014] (a homotopybased method to find the spectrum of delay equations is described in Surya et al. [2017]). By increasing the number of POD modes gives a more accurate approximation for the rightmost eigenvalues. Near the eigenvalues the magnitude of the characteristic function (21) is small. Thus substituting the computed eigenvalues into (21) should result in a residual of small magnitude. This is shown in Fig. (6) for increasing number of POD modes. 4.2 DMD results

We demonstrate the methods described in Sec. 2 and Sec. 3 on the Hayes equation (Hayes [1950]), which is one of the simplest delayed system in the following form dx(t) = ax(t) + bx(t − τ ). (20) dt The infinitely many eigenvalues of the Eq. (20) satisfy the characteristic equation: (21) λ − a − be−λτ = 0. We numerically solved the Eq. (20) with parameters a = −3, b = −3.5 and τ = 1 for t ∈ [0, 10] and initial function 255

When the pseudo-inverse X + 1 calculated, a rank reduced A ˜ (denoted by A) can be computed by truncating the matrices U 1 , V 1 and Σ1 . For example, a rank-4 approximation of A is given by   0.9971 −0.0800 −0.0215 0.0371 0.0464 0.9950 −0.0850 0.0057  ˜= A (23)  0.0004 −0.0008 0.9324 0.2509  . −0.0003 0.0005 −0.1454 0.9660 The computation of the DMD modes ψi (which come in conjugate transpose pairs) is described in Section 1.4 of

2018 IFAC TDS Budapest, Hungary, June 28-30, 2018 256

Balázs Heizer et al. / IFAC PapersOnLine 51-14 (2018) 254–258

Fig. 5. Comparison of the computed eigenvalues Fig. 2. The first four POD modes

Fig. 6. Residual at the approximated rightmost eigenvalue

Fig. 3. The ai (t) functions corresponding to the first 4 POD modes

Fig. 7. Real part of the first two DMD mode pairs Kutz et al. [2016]. The real part of the first two DMD mode pairs and the real part of their amplitudes are shown in Figs. (7) and (8), respectively. The numerical solution xnum (t) and its approximation computed by DMD are ˜ are approximating shown in Fig. (9). Eigenvalues of A the eigenvalues of Eq. (20) (Fig. (10)). By increasing ˜ the DMD gives better approximation the rank of A, for the rightmost eigenvalues. Near the eigenvalues the magnitude of the characteristic function (21) is small. Thus substituting the computed eigenvalues into (21) should result in a residual of small magnitude. This is shown in Fig. (11) for increasing number of POD modes.

Fig. 4. Numerical solution and POD approximation

256

2018 IFAC TDS Budapest, Hungary, June 28-30, 2018

Balázs Heizer et al. / IFAC PapersOnLine 51-14 (2018) 254–258

257

Fig. 11. Residual at the approximated rightmost eigenvalue

Fig. 8. The amplitudes of the first two DMD mode pairs

Fig. 12. The stability chart of the Hayes equation, colored by the number of POD modes used to reconstruct the original data less than 0.1% error. The thick dashed line is the stability boundary.

Fig. 9. Numerical solution and DMD approximation

experiments show that our findings apply well for other parameter pairs. Further investigations in the range of applicability and error estimation will be the topic of a full research paper. There is an interesting connection between the energy content of the POD modes and the stability of the Hayes equation. We found that in the stable region more POD modes were needed to reach a certain energy content. We attribute this to the fact that in the unstable region the most unstable modes overwhelm the others. The region where already one mode captured the behavior of the system, we found a single positive real leading eigenvalue, while where the dynamics was approximated by two modes, the leading eigenvalues were a complex conjugate pair. In Fig. (12) we depict how many POD modes needed to reconstruct to 99.99% of the energy content for a given (a, b) parameter pair.

Fig. 10. Comparison of the computed eigenvalues 5. DISCUSSION AND CONCLUSION Two data-driven methods were presented to compute the eigenvalues (and rightmost eigenvalues in particular) of delay-differential equation based on snapshots of the numerical solution xnum (t). This was shown by solving the Hayes equation for a parameter pair (a, b), but numerical 257

2018 IFAC TDS 258 Budapest, Hungary, June 28-30, 2018

Balázs Heizer et al. / IFAC PapersOnLine 51-14 (2018) 254–258

REFERENCES Amann, A., Sch¨ oll, E., and Just, W. (2007). Some basic remarks on eigenmode expansions of time-delay dynamics. Physica A: Statistical Mechanics and its Applications, 373, 191–202. Asl, F.M. and Ulsoy, A.G. (2003). Analysis of a system of linear delay differential equations. Journal of Dynamic Systems, Measurement, and Control, 125(2), 215–223. Bellen, A. and Maset, S. (2000). Numerical solution of constant coefficient linear delay differential equations as abstract cauchy problems. Numerische Mathematik, 84(3), 351–374. Breda, D., Maset, S., and Vermiglio, R. (2014). Stability of Linear Delay Differential Equations: A Numerical Approach with MATLAB. Springer. Brunton, B.W., Johnson, L.A., Ojemann, J.G., and Kutz, J.N. (2016). Extracting spatial–temporal coherent patterns in large-scale neural recordings using dynamic mode decomposition. Journal of neuroscience methods, 258, 1–15. Hayes, N. (1950). Roots of the transcendental equation associated with a certain difference-differential equation. Journal of the London Mathematical Society, 1(3), 226– 232. Hua, J.C., Roy, S., McCauley, J.L., and Gunaratne, G.H. (2016). Using dynamic mode decomposition to extract cyclic behavior in the stock market. Physica a: Statistical mechanics and its applications, 448, 172–180. Kutz, J.N. (2013). Data-driven modeling & scientific computation: methods for complex systems & big data. Oxford University Press. Kutz, J.N., Brunton, S.L., Brunton, B.W., and Proctor, J.L. (2016). Dynamic mode decomposition: data-driven modeling of complex systems, volume 149. SIAM. Michiels, W., Jarlebring, E., and Meerbergen, K. (2011). Krylov-based model order reduction of time-delay systems. SIAM Journal on Matrix Analysis and Applications, 32(4), 1399–1421. Proctor, J.L. and Eckhoff, P.A. (2015). Discovering dynamic patterns from infectious disease data using dynamic mode decomposition. International health, 7(2), 139–145. Schmid, P.J. (2010). Dynamic mode decomposition of numerical and experimental data. Journal of fluid mechanics, 656, 5–28. Surya, S., Vyasarayani, C., and Kalm´ ar-Nagy, T. (2017). Homotopy continuation for characteristic roots of delay differential equations using the lambert w function. Journal of Vibration and Control, DOI: 10.1077546317717629. Tu, J.H., Rowley, C.W., Luchtenburg, D.M., Brunton, S.L., and Kutz, J.N. (2013). On dynamic mode decomposition: theory and applications. arXiv preprint arXiv:1312.0041.

258