ON THE MONTE CARLO SIMULATION OF. MOMENT LYAPUNOV EXPONENTS

Wei-Chau Xie and Qinghua Huang

Department of Civil Engineering,

University of Waterloo, Canada
E-mail: xie@uwaterloo. ca, q2huang@engmail. waterloo. ca

Abstract

The moment Lyapunov exponents are important characteristic numbers for determining the dynam­ical stability of stochastic systems. Monte Carlo simulations are complement to the approximate analytical methods in the determination of the moment Lyapunov exponents. They also provide criteria on assessing how accurate the approximate analytical methods are. For stochastic dynamical systems described by Ito stochastic differential equations, the solutions are diffusion processes and their variances may increase with time of simulation. Due to the large variances of the solutions and round-off errors, bias errors in the simulation of momemt Lyapunov exponents are significant in the cases of improper numerical approaches. The improved estimation for some systems is presented in this paper.

1 Introduction

The moment Lyapunov exponents, which are define by

Л(р)= lim±logE [||X(t)f] , (1)

t

characterize the moment stability of a stochastic dynamical system with state vector X(f), where E[•] denotes expectation and ||-|| denotes a suitable vector norm. The pth moment of the response of the system is asymptotically stable if Л(р) < 0. Moreover, Л'(0) is equal to the largest Lyapunov exponent A, which is defined by

A = lim log ||X(f)|| (2)

and describes the almost-sure or sample stability of the system.

Even if the solution of a system is almost-sure stable, i. e. ||X(f) || ^0 as t ^ to w. p.l at a negative exponential rate, the chance that ||X(f)|| takes large values may lead to the instability of the pth moment. Hence, it is important to obtain the moment Lyapunov exponents such that the complete properties of dynamic stability of stochastic systems can be described.

Although moment Lyapunov exponents may be obtained by approximate analytical methods, such as stochastic averaging or perturbation, in some cases, Monte Carlo simulations have to be applied when it is difficult to do so. On the other hand, the accuracy of the approximate analytical methods needs to be verifed by numerical simulations.

A numerical algorithm for determining the moment Lyapunov exponents using Monte Carlo simula­tion developed by Xie (Xie, 2005) is described briefly below.

Consider a dynamical system whose state vector satisfies the following Ito stochastic differential equation

dX(f) = m(X, f)df + a(X, f)dW(f), (3)

627

M. Pandey et al. (eds), Advances in Engineering Structures, Mechanics & Construction, 627-636.

© 2006 Springer. Printed in the Netherlands.

where W(i) is a vector of standard Weiner processes. Define

where Tv = КД is the time interval for normalization, Д is the time step of iteration for an appropriate discrete scheme used to solve equation (3) numerically. In order to avoid data overflow or underflow, when system (3) is linear, at time rnTy the sth sample of response is normalized using

One revision to correct the insufficiency of the above algorithm is to normalize the response by their expectation. For a linear system, by defining

The solution of equation (3) is a diffusion process and its variance may increase significantly with time. Although equation (10) is exact theoretically when M is large enough, there are two main reasons which will lead to significant numerical errors. According to the Central Limit Theorem (Feller, 1965), for independent and identically distributed (i. i.d.) random variables x , • •• with the same mean p and variance a2, the sample average x = ^ Y^=i хг will ten^ to fhe normal distribution N(p, a2/n). This means that equation (7) will not give acceptable results of the expected values in the cases when the variance of response is so large that it is impossible to reduce the error of estimation with a finite sample size. On the other hand, due to the finite lengths of floating-point representations in computer, when two numbers are summed up, the smaller one will be neglected if the difference of their exponent bits exceeds the limit. Thus this truncated error in estimating the expectations will be dominant for simulations with large variances even if the chance that the responses take extremely large values is rare.

A simple example is the first-order linear stochastic system

dx(t) = ax(t)dt + ax(t)dW (t), (11)

where a and a are real constants. Its solution with the initial condition ж(0) = 1 is

x(t)= e(a-Ї )t+aW W. (12)

The moment Lyapunov exponents of this system are given by

Л(р) = | {p – l)^2 + 2a] , (13)

and the variance of norm is

Var [x(t)]= e2at(e*2f – 1). (14)

It can be seen that the variance will increase exponentially with time when a ^0.

P

Figure 1: Moment Lyapunov exponents for different times of simulation

Figure 1 shows the numerical results using the explicit Euler scheme for different times T in the case a = 0, a =1. The time step for iteration is Д = 0.001, the sample size is 5000 and equation (10) is used to determine the approximate moment Lyapunov exponents. It is obvious that the longer the time for simulation, the worse the results.