Finite size Lyapunov exponent at a saddle point

Applied Mathematical Modelling 39 (2015) 4523–4533

Contents lists available at ScienceDirect

Applied Mathematical Modelling journal homepage: www.elsevier.com/locate/apm

Finite size Lyapunov exponent at a saddle point Leonid I. Piterbarg Department of Mathematics, University of Southern California, Kaprielian Hall, Room 108, 3620 Vermont Avenue, Los Angeles, CA 90089-2532, United States

a r t i c l e

i n f o

Article history: Received 8 May 2013 Received in revised form 6 December 2014 Accepted 5 January 2015 Available online 14 January 2015 Keywords: FSLE Saddle point Hyperbolic stochastic system

a b s t r a c t A simple stochastic system is considered modeling Lagrangian motion in a vicinity of a hyperbolic stationary point in two dimensions. We address the dependence of the Finite Size Lyapunov Exponent (FSLE) k on the diffusivity D and the direction of the initial separation h. It is shown that there is an insignificant difference between the curves k ¼ kðhÞ for pure dynamics (D ¼ 0) and for infinitely large noise (D ¼ 1). For small D a well known boundary layer asymptotic is employed and compared with numerical results. Ó 2015 Elsevier Inc. All rights reserved.

1. Introduction In recent years FSLE has become a popular tool for investigating mixing in ocean and atmospheric flows, e.g. [1,2]. In theoretical works the focus was mostly on the scaling of kðdÞ as a function of the initial separation magnitude d for flows close to isotropic, e.g. [3–6]. In the presented work we consider an extremely anisotropic flow and first suggest a closed solution for the case of infinitely large noise. It is found that the anisotropy in k is almost not affected by an infinite diffusivity. To analyze the case of small D we use the boundary layer approach, e.g. [7], as well as some heuristic arguments, and compare results with direct numerical computations of k. Next we proceed to exact formulations. Let Xðt; aÞ be a Lagrangian trajectory starting at a in a stochastic velocity field uðt; xÞ

X_ ¼ uðt; XÞ;

Xð0Þ ¼ a;

ð1Þ

and define the separation process

Zðt; a; DaÞ ¼ Xðt; a þ DaÞ  Xðt; aÞ: Let

s ¼ sða; Da; aÞ be the first moment when the separation hits ajDaj where a > 1 is a prescribed threshold

s ¼ infft > 0 : jZðtÞj P ajDajg: Following to [3,4] define the forward FSLE by

kþ ¼

ln a ; h si

where the angels mean the ensemble averaging, and in a similar way define the backward FSLE k by replacing in (1) uðt; xÞ with uðt; xÞ. Finally, FSLE itself is defined by

k ¼ kða; Da; aÞ ¼ kþ  k : http://dx.doi.org/10.1016/j.apm.2015.01.002 0307-904X/Ó 2015 Elsevier Inc. All rights reserved.

4524

L.I. Piterbarg / Applied Mathematical Modelling 39 (2015) 4523–4533

The problem of analytical studying k for a general stochastic velocity field is extremely hard, in particular because the separation equation

Z_ ¼ uðt; Xðt; a þ DaÞÞ  uðt; Xðt; aÞÞ; is not closed. To advance we make a few assumptions under which it allows for a closure. First, assume

uðt; xÞ ¼ UðxÞ þ u0 ðt; xÞ;

ð2Þ

where UðxÞ ¼ Gx is a deterministic velocity field with constant matrix G (shear tensor) and u0 ðt; xÞ is a stochastic field with zero mean. Conditions for representation (2) can be found for example in [8]. From (2) it follows that

Z_ ¼ GZ þ u0 ðt; Xðt; a þ DaÞÞ  u0 ðt; Xðt; aÞÞ: Next assume that the stochastic (turbulent) component u0 ðt; xÞ is homogeneous in space and u0 ðt; x1 Þ; u0 ðt; x2 Þ are statistically independent whenever x1 – x2 , i.e. it is a white noise in space. These assumptions imply that statistical characteristics of the Lagrangian separation velocity

v0 ðtÞ ¼ u0 ðt; Xðt; a þ DaÞÞ  u0 ðt; Xðt; aÞÞ; depend on time only and we arrive at a Langevin equation

Z_ ¼ GZ þ v 0 ðtÞ:

ð3Þ

The main goal of this work is to investigate k for the simplest hyperbolic structure of the shear tensor

G¼



l 0 l

0

 ;

implying flow incompressibility, under additional assumption that v 0 ðtÞ is a Gaussian white noise with uncorrelated components. In other words we address FSLE for the following Eulerian velocity field

uðt; x; yÞ ¼ lx þ

pffiffiffiffi Dw_ 1 ;

v ðt; x; yÞ ¼ ly þ

pffiffiffiffi Dw_ 2 ;

ð4Þ

with independent white noises, where D is the diffusivity and l the Lyapunov exponent for the corresponding hyperbolic dynamical system. Regarding the case of large diffusivity we address an arbitrary constant shear tensor. Originally, the problem of investigating FSLE was motivated by physical oceanography applications, in particular a question of interest was how effective is FSLE in detecting hyperbolic and other stagnation points of the large scale circulation in the presence of intense small scale turbulence. In view of the strong assumptions leading to (4) that model is hardly adequate to real ocean turbulence in the upper layer. However, it is presumably the only exactly solvable model reflecting both hyperbolic circulation and stochastic perturbations. So, all the effects derived from (4) would be of some interest for further studies in frameworks of more realistic models. First we briefly discuss the case D ¼ 0 for (4), then investigate large and small D, and finally summarize conclusions. 2. Pure dynamics (D ¼ 0) In this case an equation for the separation time

x20 e2ls

þ

y20 e2ls

2

¼a

ðx20

þ

s reads

y20 Þ;

where ðx0 ; y0 Þ is the initial separation. Substituting x0 ¼ r0 cos h; y0 ¼ r0 sin h and solving for

s get

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi! 1 a2  a4  sin2 2h : s¼ ln 2l 2 cos2 h Sign þ should be taken otherwise s < 0. The separation time for the opposite velocity is obtained by replacing h with h þ p=2 Hence

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi! pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi! 1 a2 þ a4  sin2 2h 1 a2 þ a4  sin2 2h  T ¼ ln ln T ¼ ; 2 2l 2l 2 cos2 h 2 sin h þ

and FSLE for the unperturbed system

k0 ¼ ln a



 1 1   ; Tþ T

ð5Þ

does not depend on the reference point a nor on the magnitude of initial separation jDaj because of translation invariance. Thus it is completely defined by a and the angle h between Da and the x-axis. Notice that for horizontal initial separation

L.I. Piterbarg / Applied Mathematical Modelling 39 (2015) 4523–4533

4525

(h ¼ 0) corresponding to the largest rate of divergence we have k ¼ l and for vertical separation (h ¼ p=2), k ¼ l. Another important observation is that for p=4 < h 6 p=2 the limit of T þ as a ! 1 is not zero since for such h’s the separation first decreases during time

t0 ¼

ln tan h : 2l

A special point of interest is the limit of k as a ! 1. It can be found that

limk ¼ l cos 2h; a!1

which follows from the expression for the derivative of the forward separation time at a ¼ 1 þ

T a ja¼1 ¼

8 1 > < cos 2h > :

h < p4 ;

 cos12h

hP

p 4

þ

and the mentioned fact that T remains finite for

p=4 < h 6 p=2 as a ! 1

3. Separation time for model (4) and FSLE for large D In general case D > 0 the separation Eq. (3) for system (4) in coordinate form becomes

X_ ¼ lX þ

pffiffiffiffiffiffiffi 2Dw_ 1 ;

Y_ ¼ lY þ

pffiffiffiffiffiffiffi 2Dw_ 2 ;

where Z ¼ ðX; YÞ. The process ðX; YÞ is a Markov process with the generator

L ¼ DD þ lx@ x  ly@ y ; where D is Laplacian and @ x ; @ y are the partial derivatives in x and y respectively. It is well known, e.g. [9], that under wide conditions the mean exit time T ¼ Tðx; yÞ from a bounded region @G satisfies the following zero boundary condition problem

LT ¼ 1;

Tj@G ¼ 0;

where @G is the boundary of G assumed to be smooth and ðx; yÞ 2 G is the starting point of a particle described by the Markov process with generator L. Let d > 0 be a certain space scale and Tðx; yÞ ¼ Tðx; y; adÞ be the mean time for the separation to hit ad with the initial pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi separation ðx; yÞ; r < ad, where r ¼ x2 þ y2 . The above general theory leads to

DDT þ lxT x  lyT y ¼ 1;

Tjr¼ad ¼ 0:

With that notation

kþ ¼

ln a ; Tðd cos h; d sin h; adÞ

where d is the magnitude of initial separation and h the polar angle. Thus, the separation time s still does not depend on the referenced point because of the translation invariance of (4). However with the noise k depends on both the magnitude and direction of the initial separation as well as on a. After proceeding to dimensionless variables

x ! x=d;

y ! y=d;

T ! lT;

we arrive at

D0 DT þ xT x  yT y ¼ 1;

Tjr¼a ¼ 0;

ð6Þ

where 0

D ¼

sffiffiffiffiffiffiffiffi D d2 l

;

is the dimensionless diffusivity which henceforth will be again denoted by D In polar coordinates the solution of (6) is expressible in the form (see Appendix)

Tðr; hÞ ¼ e

r2 cos 2h 4D

 2 1 X r cosð2nhÞ; gn 4D 1

ð7Þ

4526

L.I. Piterbarg / Applied Mathematical Modelling 39 (2015) 4523–4533

where

g n ðzÞ ¼ K n ðzÞ

Z

z

0

I2n ðyÞdy  In ðzÞ

Z

z

In ðyÞK n ðyÞdy þ CðsÞIn ðzÞ:

0

K n ; In are the modified Bessel functions, s ¼ a2 =4D, and CðsÞ is a constant ensuring the boundary condition. An exact expression for C is given in Appendix. Let D ! 1, then s ! 0 and z ¼ 1=4D ! 0. It can be shown that 2

g 0 ðzÞ ¼ zða 21Þ þ Oðz3 ln zÞ g 1 ðzÞ ¼ z

2 ða2 1Þ

g n ðzÞ ¼ O

3



znþ1 2n n!

þ Oðz3 Þ



;

n P 2:

Thus for large D the series (7) uniformly converges and one gets for T þ ¼ Tð1; hÞ

Tþ ¼

  zða2  1Þ 2 1  z cos 2h þ Oðz2 ln zÞ ; 2 3

and similarly for the adjoint problem

T ¼

  zða2  1Þ 2 1 þ z cos 2h þ Oðz2 ln zÞ : 2 3

Despite both kþ and k indefinitely increase as D ! 1 the following limit

k1 ¼ lim k ¼ D!1

8 ln a cos 2h; 3ða2  1Þ

is finite. Thus, an important conclusion is that in the limiting case of large noise the sharp anisotropic structure of k is conserved and in fact the dependence on h has the same tendency as in the case D ¼ 0 as illustrated in Fig.1. Notice that symmetries of the system (4) imply

T þ; ðhÞ ¼ T þ; ðp  hÞ;

T þ ðp=2  hÞ ¼ T  ðhÞ;

thereby k is symmetric about p=4. To sum up let us go back to dimensional variables. The separation time is given by

2

D



2d4 l cos 2h 3D

!

2

  þ O D3 :

1

α=1, D=0 α=1, D=∞ α=2, D=0 α=2, D=∞ α=5, D=0 α=5, D=∞

0.8 0.6 0.4 0.2

λ

T¼

a2  1 d2

0

−0.2 −0.4 −0.6 −0.8 −1 0

0.5

θ

1

1.5

Fig. 1. k vs h for different a at zero diffusivity (curves with no markers) and infinite diffusivity (markers).

L.I. Piterbarg / Applied Mathematical Modelling 39 (2015) 4523–4533

4527

So, in the first order in D1 the separation time is indeed determined by the diffusivity only while the dynamics is merely manifested in the second order. However, just this term survives when one switches to k because the first term is the same for kþ and k

k1 ¼

4 ln a l cos 2h; 3ða2  1Þ

and that limit is completely determined by the dynamics. In fact, the large diffusivity limit can be found for any stagnation point of a linear system. Indeed, let us represent the solution of the Eq. (6) in which lx and ly are replaced with

v ðt; x; yÞ ¼ cx þ dy;

uðt; x; yÞ ¼ ax þ by; respectively, as

T¼

1 0 1 T þ 2 T1; D D

and hence

k1 ¼ ln a lim

D!1

1 1 0 T D

þ D12 T

1 1 D

1 T 0  D12 T 1

!

2T 1 ln a ¼   2 ; T0

T 0 and T 1 are the solutions of

DT 0 ¼ 1;

DT 1 þ uT 0x þ v T 0y ¼ 0;

corresponding to the zero boundary conditions at r ¼ a. Not difficult to find

T0 ¼

a2  r2 4

;

T1 ¼

ða þ dÞðr 4  a4 Þ r 2 ða  dÞðr 2  a2 Þ r2 ðb þ cÞðr 2  a2 Þ cos 2h þ sin 2h; þ 32 48 48

and hence

k1 ¼ ca

  ða þ dÞð1 þ a2 Þ 2ða  dÞ 2ðb þ cÞ cos 2h þ sin 2h ; þ 2 3 3

ca ¼

ln a : a2  1

Of course, the perturbation approach is always limited, however because it gives the same result for the saddle point as the rigorous derivation does, one can believe that the last expression is true in general. Expressions for k1 =ca are shown in Fig.2. for some particular cases: Gyre for which ða; b; c:dÞ ¼ ð0; 1; 1; 0Þ, Saddle (1; 0; 0; 1), Shear (0; 0; 1:0), Divergence point (1; 0; 0; 1), and Convergence point (1; 0; 0; 1). 4. Small D Series (7) is of little help for studying small diffusivity since its terms go to infinity. Instead we return to the basic Eq. (6) written in polar coordinates and change notation D ¼ e to match the goal

eDT þ LT ¼ 1; TjC ¼ 0; where L ¼ x@ x  y@ y ¼ r cos 2h@ r  sin 2h@ h , C ¼ fðr; hÞ : G ¼ fðr; hÞ : r < a; 0 < h < p=2g.

ð8Þ r ¼ a; 0 < h < p=2g is the boundary of the region of interest

Fig. 2. Normalized FSLE at infinite diffusivity k1 =ca for different types of circulation around a stagnation point.

4528

L.I. Piterbarg / Applied Mathematical Modelling 39 (2015) 4523–4533

First, notice that the solution of (8) converges to the separation time

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi! 1 a2 þ a4  r4 sin2 2h T ¼ ln ; 2 2r 2 cos2 h 0

ð9Þ

for the unperturbed system as e ! 0 [10, Section 2.3]. Then, applying the method of boundary layer decomposition [7] one can obtain the first order asymptotic

~ hÞ þ Oðe2 Þ; Tðr; hÞ ¼ Tðr; of the solution of (8) in the form

   ~ hÞ ¼ T 0 ðr; hÞ þ eT 1 ðr; hÞ  exp  S T 0 ða; hÞ þ eT 1 ða; hÞ þ b1 S þ 1 b2 S2 ; Tðr;

e

e

ð10Þ

where T 0 ðr; hÞ is given by (9) and it satisfies LT 0 ¼ 1 with zero boundary conditions only on the part

C1 ¼ fðr; hÞ : r ¼ a; 0 < h < p=4g of C. The first order regular correction T 1 satisfies LT 1 ¼ DT 0 ;

ð11Þ

with zero boundary conditions on C1 as well and is given by

T 1 ¼ Bðr; hÞ  Bða; bðr; hÞÞ; 1

with b ¼ 12 sin



r 2 sin 2h



a2

ð12Þ

and

Bðr; hÞ ¼ r2 Aðr 2 sin 2hÞ cos 2h þ

1 ; 2r 2 cos2 h

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

AðtÞ ¼

a6 þ ða4 þ t2 Þ a4  t2

: pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ða2 a4  t 2 þ a4  t 2 Þ a4  t2

Representation of T 1 in form (12) allows for direct verification of (11) and the boundary conditions accounting for

DT 0 ¼ 2r 2 Aðr 2 sin 2hÞ þ

1 ; r 2 cos2 h

and an obvious observation that Lf ¼ 0 for any smooth f ¼ f ðr 2 sin 2hÞ. The exponent is given by

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi! 1 a2 þ a4  r4 sin2 2h  4 4 2  S ¼ 2 ln a  r sin 2h : 2 2a 2r 2 sin h The meaning of the first factor is the time required for a particle driven by the unperturbed system to reach ðr; hÞ starting at the point ða; hÞ on the second part C2 ¼ fðr; hÞ : r ¼ a; p=4 < h < p=2g of the boundary C, thereby S ¼ 0 on C2 . The second factor is coming from Laplacian in (8). Notice the boundary value for h > p=4

T 1 ða; hÞ ¼

8 cos4 h  8 cos2 h þ 1 2

2a2 sin h cos2 h cos2 2h

;

which is singular at h ¼ p=4; p=2. Finally, explicit expression for prefactors b1 ðr; hÞ; b2 ðr; hÞ from (10) can be found in Appendix. All in all, the meaning of approximation (10) is that it satisfies, first, Eq. (8) up to the order of e2 , second, zero boundary conditions on C2 exactly, and, finally, zero boundary conditions on C2 approximately with an exponentially small error of order expða2 cos2 2h ln cot h=eÞ. It follows from the last expression that the approximation may not be satisfactory in some vicinity of h ¼ p=4 which indeed was observed in numerical experiments when a close to 1. Details of derivation of (10) are given in Appendix. In first two panels of Fig. 3 we compare the separation time for unperturbed system T 0 ð1; hÞ (blue) and the boundary layer ~ hÞ (black) with the exact solution Tð1; hÞ (red) obtained from series (7) for particular values D ¼ 0:05; a ¼ 2 in asymptotic Tð1; the range 0 < h < 7p=16. A validation of the computations via the series (7) is given in Appendix. As one can see the asymptotic steadily improves the zero approximation T 0 up to h  3p=8, after which both approximations become equally unsatisfactory. On the interval ð7p=16; p=2) they both make little sense because of blowing up at h ¼ p=2. Moreover T~ turns out to be even worse than T 0 ð1; hÞ. Notice that for a well above 1 the difference between boundary layer approximation T~ and the regular approximation, given by the first two terms in (10), is almost invisible. However, as a approaches 1 the picture becomes more sophisticated: T~ better approaches T only outside some vicinity of h ¼ p=4 while for other h it is worse than the regular approximation.

4529

L.I. Piterbarg / Applied Mathematical Modelling 39 (2015) 4523–4533 D=0.05,α=2

D=0.05,α=2

2.8

0.5

3 Numerical Unperturbed Asymptotic

2.5

Asymptotic,α=1.9

Unperturbed Asymptotic

0.4

Numerical,α=1.5

2.5

0.2

1.5

Asymptotic,α=1.5

2.6

0.3 2

Numerical,α=1.9

2.7

0.1

2.4

0

2.3

−0.1

2.2

−0.2

2.1

1

0.5 0

0.2

0.4

0.6

0.8

θ

1

−0.3 0

1.2

0.2

0.4

0.6

θ

0.8

1

2 0.05 0.055 0.06 0.065 0.07 0.075 0.08 0.085 0.09 0.095

1.2

0.1

D

Fig. 3. Small diffusivity. Left: T (numerical from (7)), T 0 (unperturbed given by (9)), and T~ (asymptotic (10)) vs. h on interval ½0; 7p=16 for D ¼ 0:05; a ¼ 2 Center: Approximation errors T 0  T; T~  T (right) for the same parameters. Right: Comparison Tð1; p=2Þ with one-dimensional asymptotic (15) for different a and D.

Since (10) fails at h ¼ p=2 we tried an alternative approach exploiting a one dimensional approximation. Namely, near h ¼ p=2 it seems reasonable to neglect the dependence on y when a is essentially greater than 1 and consider the following one-dimensional problem

^ eT^ xx þ xT^ x ¼ 1; Tj x¼a ¼ 0:

ð13Þ

^ The solution TðxÞ of (13) is given by

^ TðxÞ ¼2

Z

pffiffiffiffi

a= 2e

pffiffiffiffi x= 2e

uðyÞUðyÞdy;

ð14Þ

where 2

uðyÞ ¼ ey ; UðyÞ ¼

Z

y

2

ez dz:

0

Then we suggest the following approximation

 p a ^ ¼ ln pffiffiffiffiffiffi þ Oð1Þ: T 1;  Tð0Þ 2 2e

ð15Þ

^ (black) in the third panel of Fig.3 for modestly small D ¼ e and two values of a well We compare Tð1; p=2Þ (red) and Tð0Þ separated from 1: a ¼ 1:9 (no marks) and a ¼ 1:5 (with marks. Thus, (15) does really work for relatively large a while for a close to 1 it turns out unsatisfactory. Regarding to FSLE we first compared both k0 given by (5) and the boundary layer approximation ~ k based on (10), with numerical k obtained from (7). In first two panels of Fig.4 such a comparison is shown for a ¼ 2; D ¼ 0:05 in the whole range

D=0.05, α=2

D=0.05, α=2

1

1

0.3 Numerical Unperturbed Asymptotic

0.8 0.6

α=1.1

Unperturbed Asymptotic

α=1

0.2

0.998

2

1−4D

1−0.16D

0.997

0.4

0.1 λ

0

0

1−0.38D

0.996

0.2

0.995 0.994

−0.2 −0.1

−0.4 −0.6

0.993 0.992

−0.2

−0.8 −1

α=1.5

0.999

0.991

0

0.5

1

θ

1.5

0

0.5

1

θ

1.5

0.99 0.002 0.004 0.006 0.008 0.01 0.012 0.014 0.016 0.018 0.02

D

~ (black) vs h for a ¼ 2 Fig. 4. Small diffusivity. Left: k obtained from (7) (red), FSLE for unperturbed system k0 (blue), and boundary layer approximation k and D ¼ 0:05. Center: The error of approximation of k by k0 (blue) and ~ k (black) for a ¼ 2 and D ¼ 0:05. Right: Approximation of k at h ¼ 0 via (16) vs D for different a. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.).

4530

L.I. Piterbarg / Applied Mathematical Modelling 39 (2015) 4523–4533 0.35

0.8

θ=π/8 θ=π/4 θ=π/2

0.3

Simulation Numerical

0.7 0.6

0.25

0.5 λ

0.2

0.4

0.15

0.3 0.1

0.2

0.05 0

0.1

50

52

54

56

58

60

N

0 0

0.1

0.2

0.3

0.4 θ

0.5

0.6

0.7

Fig. 5. The error of (7) vs cut off N for D ¼ 0:001; a ¼ 1:1 (left). Comparison of k obtained by simulations and numerics with N ¼ 10 for D ¼ 1; a ¼ 1:1 (right).

k is even worse than k0 because of ð0; p=2Þ. As one can see both approximations are not satisfactory near the endpoints while ~ the mentioned trouble with h ¼ p=2 for the separation time asymptotics. However in the range ð0:16p=2; 0:84p=2Þ the former is better than the latter even though not significantly (central panel). Finally, we zoom in k at a key value h ¼ 0. For that let us neglect k at this point which is basically the reciprocal of T þ ð1; p=2Þ and use only the regular part of approximation (12) for T þ . The result reads

k  1  ch2  c1 e  c2 eh2 ;

c¼

2 a4  1 ða2  1Þ ða2  1Þð4a2  1Þ ; c1 ¼ ; c2 ¼ : 4 4 2a ln a 4a6 ln a 2a ln a

In particular for h ¼ 0

k  1  c1 ðaÞe:

ð16Þ

In the right panel we compare (16) with the direct computations via (7) a ¼ 1:1 and a ¼ 1:5 with c1 ¼ 0:16 and c1 ¼ 0:38 respectively. Here we also show kðDÞ for a ¼ 1 when c1 ¼ 0 and it should be compared with the second order approximation k ¼ 1  4e2 (red curve). Computationally the case a ¼ 1 was investigated by using the derivative of (7) in a cos 2h 4e

T a ð1; hÞja¼1 ¼ e

1 X 1

1 In ð1=4eÞ

Z 0

1=4e

I2n ðyÞdy cosð2nhÞ:

ð17Þ

Finally, let us compare (16) with the one-dimensional approximation (14). Namely, the latter implies

^k ¼

ln a ffi ; R a=pffiffiffi 2e y2 R y z2 p ffiffiffi ffi 2 e e dzdy 1= 2e

and for small

0

e

k  kþ  ^k  1  c1 ðaÞe; c1 ðaÞ ¼

a2  1 : 2a2 ln a

Since c1 ð1Þ – 0 we conclude that the 1D approximation is not able to adequately describe the behavior of k at h ¼ 0 in contrast to h ¼ p=2. 5. Conclusions A simplest linear hyperbolic system perturbed by spatially uncorrelated white noise was considered. An explicit solution was found for a partial differential equation covering the mean separation time of two Lagrangian trajectories. That solution was used to find the limit of FSLE as diffusivity indefinitely increases. This limit turns out to be finite and it well reproduces the FSLE dependence on the direction h of the initial separation for the unperturbed system. The same limit was computed for a general linear velocity field by perturbation means. To investigate small diffusivity asymptotics we first applied a general result from [10] to establish that the limit of FSLE as D approaches zero coincides with that of the unperturbed system. Then, we made use of the boundary layer asymptotic for elliptic problems [7]. Comparing that asymptotic with numerical results demonstrate its helpfulness in a certain range of the polar angle separated from the endpoints h ¼ 0 and h ¼ p=2.

L.I. Piterbarg / Applied Mathematical Modelling 39 (2015) 4523–4533

4531

The suggested asymptotics are not satisfactory for the threshold a close to 1 because of a singularities at h ¼ p=4; p=2 of the boundary layer approximation. Thus, investigating the behavior of FSLE when both D and a  1 approach zero is left for the future. Another important direction of further research would be extending the results to a more general shear tensor in the Langevin Eq. (3). We restricted ourselves mostly to the diagonal tensor because only in this case the variables r; h in (6) can be separated implying the exact expansions (7, 17) needed for testing the suggested asymptotics. In summary, our results show that FSLE is an extremely efficient instrument for detecting saddle points of dynamical systems regardless of the intensity of the stochastic perturbations. Acknowledgments The support of the Office of Naval Research under Grant N00014–11-1–0369 and NSF under Grant CMG-1025453 is greatly appreciated. The author thanks Sergey Lototsky, Lenya Ryzhik and Alexei Novikov for fruitful discussions. Comments of anonymous reviewers helped to improve presentation. Appendix A A.1. Derivation of (7) Substituting

Tðx; yÞ ¼ e

y2 x2 4D

f ðx; yÞ;

into (6) get

Df 

x2 þ y 2

f ¼

2

4D

1 x2 y2 e 4D : D

This equation allows variable separation in polar coordinates accounting for [11, 9.6.34] 1 X

ea cos 2h ¼

In ðaÞ cosð2nhÞ:

1

Namely 1 X

f ðr; hÞ ¼

f n ðrÞ cosð2nhÞ;

1

where 00

0

r2 f n þ rf n 



  2 r2 r 2 f ; þ 4n ¼  I n n D 4D 4D2 r4

f n ðaÞ ¼ 0:

ð18Þ

Variable change z ¼ r2 =4D reduces the Eq. (18) to a standard Bessel equation for gðzÞ ¼ f ðrÞ

  z2 g 00n þ zg 0n  z2 þ n2 g n ¼ zIn ðzÞ;

g n ðaÞ ¼ 0;

a¼

a2 4D

;

with the solution bounded at z ¼ 0 given by

g n ðzÞ ¼ K n ðzÞ

Z 0

z

I2n ðyÞdy  In ðzÞ

Z

z

In ðyÞK n ðyÞdy þ CðaÞIn ðzÞ;

0

where

CðaÞ ¼

Z

a

In ðyÞK n ðyÞdy 

0

K n ðaÞ In ðaÞ

Z 0

a

I2n ðyÞdy:

A.2. Perturbation solution for (8) Following [7] introduce new coordinates ðu; tÞ by

x ¼ aet cos u;

y ¼ aet sin u:

ð19Þ

In other words t is the time for a particle to travel from a boundary point with polar coordinates ða; uÞ to the point ðx; yÞ. The inversion of (19) written in polar coordinates is given by

4532

L.I. Piterbarg / Applied Mathematical Modelling 39 (2015) 4523–4533

p 1

u ¼  sin 2

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi! 1 a2 þ a4  r4 sin2 2h t ¼ ln : 2 2 2r 2 sin h

 2  r sin 2h ; 2

1

a

2

ð20Þ

Thus G ¼ fðr; hÞj 0 < r < a; 0 < h < p=2g is mapped to fðu; tÞj p=4 < u < p=2; t > 0g and the boundaries C1 and C2 are mapped to C 1 ¼ fðu; tÞj p=4 < u < p=2; t ¼ qðuÞg and C 2 ¼ fðu; tÞj p=4 < u < p=2; t ¼ 0g respectively with qðuÞ ¼ ln tan u. In new coordinates f ðu; tÞ ¼ Tðx; yÞ satisfies

eKf þ f t ¼ 1; f jC1 [C2 ¼ 0;

ð21Þ

where

K ¼ ðt 2x þ t 2y Þ@ tt þ 2ðt x ux þ t y uy Þ@ tu þ ðu2x þ u2y Þ@ uu þ ðt xx þ t yy Þ@ t þ ðuxx þ uyy Þ@ u ; with

ux ¼

sin u t cos u t cos u t sin u t e ; uy ¼ e ; tx ¼ e ; ty ¼ e: a cos 2u a cos 2u a cos 2u a cos 2u

Dividing both parts of (21) by t2x þ t2y we lands in the conditions of Theorem 1 from [7] according to which the first order asymptotics is given by 0

1

f ðu; tÞ ¼ f ðu; tÞ þ ef ðu; tÞ þ w0 ðu; t=eÞ þ ew1 ðu; t=eÞ; where the regular perturbation is represented by the first two terms satisfying 0

0

f t ¼ 1; f jC 1 ¼ 0;

1

f t ¼ Kf

0

1

f jC 1 ¼ 0;

ð22Þ

and the last two terms come from the boundary layer expansion and satisfy

w0tt þ aðuÞw0t ¼ 0;

0

w0 jt¼0 ¼ f ðu; 0Þ; w0 jt¼1 ¼ 0;

w1tt þ aðuÞw1t þ bðuÞw0tu þ ðcðuÞt þ dðuÞÞw0t ¼ 0;

ð23Þ 1

w1 jt¼0 ¼ f ðu; 0Þ; w1 jt¼1 ¼ 0;

where

 1   a ¼ t 2x þ t2y ¼ a2 cos2 2u; b ¼ 2a t x ux þ t y uy ¼ 2 sin 2u t¼0 t¼0 :  1   2 2u 2 c ¼ @t@ t2x þ t2y jt¼0 ¼ 2a2 cos3 2u; d ¼ a txx þ t yy t¼0 ¼ 1þsin  cos 2 u cos 2u Problems (22)–(24) are easy to solve with the following results 0

f ¼ t þ qðuÞ; 1

f ¼

    1  2 Aðcos2 uÞ e2t  cot2 u  Aðsin uÞ e2t  tan2 u ; 2 2a

where

AðuÞ ¼

ðu  1Þð4u2  6u þ 1Þ uð2u  1Þ3

:

Then

w0 ¼ qðuÞeat ; and finally

  1 w1 ¼ f ðu; 0Þ  BðuÞt  CðuÞt 2 eat ; where 1

4

2

u8 cos uþ1 f ðu; 0Þ ¼ 2a82 cos ; sin2 u cos2 u cos2 2u

p B ¼ ap2 þ qa ; C ¼ 2a

: p ¼ aqðc þ 2a2 b sin 4uÞ ¼ 2a4 q cos2 2uðcos3 2u þ 2 sin 2u sin 4uÞ;   2 2 2 0 2 3 q ¼ bð2a sin 4u þ aq Þ þ aqd ¼ a cos 2u 4 sin 2u  4 cos 2u þ q cos 2uð1 þ sin 2u  cos 2uÞ Using (20) one can recompute B; C in prefactors b1 ; b2 appearing in (10). Notice that

w0t ¼ aeat ;

w0tu ¼ 2a2 ð1  atÞeat sin 4u;

ð24Þ

L.I. Piterbarg / Applied Mathematical Modelling 39 (2015) 4523–4533

4533

and hence (24) is written as

w1tt þ aw1t ¼ ðpt þ qÞeat : A.3. Validation of computations via (7, 17) To validate the computational procedure notice that for large n and fixed z, [12]

1  ez n In ðzÞ  pffiffiffiffiffiffiffiffiffi ; 2pn 2n

K n ðzÞ 

pffiffiffiffi 

p ez n : 2n 2n

Thus for a > 1

  ez n g n ðzÞ ¼ O n2 ; 2n and hence series (7) converges for each fixed z with an extremely high rate. To roughly estimate the residual of series (7) we compare it with the residual of the corresponding geometric series. Namely for N > p ¼ e=8D

RN ðzÞ ¼

X

jg n ðzÞj < 2

jnjPNþ1

1 X

pn ¼

n¼Nþ1

2pNþ1 : 1p

Thus, to ensure 99% accuracy one should cut off the series at N satisfying

2pNþ1 < 0:01SN ðzÞ; 1p where SN ðzÞ is the partial sum in (7). For small D that means simply N > e=8D, i.e. for say D ¼ 0:001 one should take N ¼ 340 which appears to be significantly overestimated. Indeed, practically it turns out that N ¼ 60 is enough as illustrated in the left panel of Fig.5. The vertical axis show the error of (7) defined as ratio jg Nþ1 ðzÞj=SN ðzÞ for three values of h. To verify (7) we also used Monte Carlo simulations averaging the separation time over 50000 samples assuming D ¼ 1 (right panel of Fig.5. Notice that series (17) converges slightly slower than for a > 1, namely its terms are of order

  ez n O n3=2 ; 2n but still fast enough to make computations efficiently. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12]

B. Joseph, B. Legras, J. Atm. Sci. 59 (2002) 1198. A. Haza, T. Ozgokmen, A. Griffa, Z. Garaffo, L. Piterbarg, Ocean Model. 42 (2012) 31–49. V. Artale, G. Boffetta, A. Celani, M. Cencini, A. Vulpiani, Phys. Fluids A 9 (1997) 3162. E. Aurell, G. Boffetta, A. Crisanti, G. Paladin, A. Vulpiani, Phys. Rev. E 53 (1996) 2337. R. Parashar, J.H. Cushman, Phys. Rev. E 76 (2007) 017201. L.I. Piterbarg, Appl. Math. Model. 36 (2012) 3464. C.J. Holland, Singular perturbations in elliptic boundary value problems, J. Differ. Equ. 20 (1976) 248–265. A.S. Monin, A.M. Yaglom, Statistical Fluid Mechanics: Mechanics of Turbulence, MIT Press, Cambridge, MA, 1975. Z. Schuss, Theory and Application of Stochastic Differential Equations, Wiley, New York, 1980. M.I. Freidlin, A.D. Wentsel, Random Perturbations of Dynamical Systems, Springer-Verlag, 1984. M. Abramovitz, I. Stegun, Handbook of Mathematical Functions, Dover, New York, 1972. Digital Library of Mathematical Functions, .