## Abstract

We present derivative-free weak and strong solutions of stochastically driven nonlinear oscillators of engineering interest using higher order forms of the locally transversal linearization (LTL) method. Unlike strong solutions, weak stochastic solutions attempt to predict only mathematical expectations of functions of the true solution and are obtainable with much less computational effort. The linearized equations corresponding to the higher order implicit LTL schemes are arrived at using backward Euler–Maruyama and Newmark expansions in order to conditionally replace nonlinear drift and multiplicative diffusion vector fields. We also briefly describe alterations through which explicit forms of such higher order linearizations are obtained. In the weak forms, the Gaussian stochastic integrals, appearing in the linearized solutions, are replaced by random variables with simpler and discrete probability distributions. The resulting local approximations to the true solution are of significantly higher formal order of accuracy, as reflected through local error estimates. Numerical illustrations are presently provided for the Duffing and Van der Pol oscillators driven by additive and multiplicative noises, which are indicative of the numerical accuracy, computational speed and algorithmic simplicity.

## 1. Introduction

Stochastically driven nonlinear oscillators, modelled using stochastic differential equations (SDEs), are of considerable importance owing to their ability to model engineering dynamical systems with uncertain inputs and parameters. Other than time-variant reliability and fatigue analyses, solutions of such SDEs are also of importance in state or parameter estimation using particle filters. When we approximate a sample path of the solution process evolving from a specific initial condition, the approximation is referred to as ‘strong’. On the other hand, a weak approximation is focused on approximating the moment of a smooth functional of the solution. In reliability analyses and filtering problems, weak approximations are preferred as a precisely pathwise nature of errors in modelling or measurements are rarely available. Additionally, one does encounter cases where weak solutions exist for SDEs while strong solutions do not. These observations emphasize the need to develop higher order weak schemes for SDEs and, in particular, noise-driven nonlinear oscillators of engineering interest. Since a more complete representation of noise is usually not available, modelling through (filtered) Wiener processes is often considered acceptable. While the Wiener process is continuous everywhere, it is nowhere differentiable (even though its formal derivative exists as a measure and is referred to as white noise), and thus complications arise while generalizing integration schemes for ordinary differential equations to SDEs. Moreover, the Wiener process has quadratic variation that is linear in time. Thus, a serious reduction of the order of accuracy occurs while developing strong or weak schemes for SDEs. However, with a weak scheme, one has the advantage of an increase in global order by *O*(*h*^{0.5}) (*h* denoting the step size) vis-à-vis its strong counterpart if the global order of the strong scheme is fractional. One way of deriving higher order schemes is to incorporate more terms (hierarchically) in the stochastic Taylor expansion; however, this would be at the expense of (possibly inaccurate) evaluations of more partial derivatives of drift and diffusion fields. Jimenez (2002) and Carbonell *et al*. (2005) have studied gradient-based local linearization schemes and their convergence for random differential equations and SDEs with Hölder continuous vector fields. However, such methods would not work for SDEs with non-differentiable fields. This highlights the need for derivative-free schemes such as stochastic Runge–Kutta (SRK) or Newmark schemes (Roy & Dash 2002, 2005*b*). An early form of the SRK scheme is by Rumelin (1982) with a strong order of 0.5 (1.0 for a rather restrictive class of SDEs). Burrage *et al*. (2004) have worked on strong explicit (Burrage & Burrage 1996, 2000) and implicit (Burrage & Tian 2004) SRK schemes. Tocino & Vigo-Aguiar (2001), Vigo-Aguiar & Tocino (2001) and Komori (2007) have also worked on various forms of SRK schemes. In the process, a maximum strong global order of 2.0 has been achieved. Roy (2001, 2004) has proposed a locally transversal linearization (LTL) for stochastically driven nonlinear oscillators and this provides a derivative-free, semi-analytical alternative to purely numerical approaches. More recently, a significant reduction in the computational effort has been achieved through a weak version of the stochastic locally transversal linearization (SLTL; Roy *et al*. submitted) with a weak global order of 1.0.

The key idea behind the LTL method is to construct a set of conditionally linear and integrable (in closed-form) vector fields over a time-step, such that with known initial conditions at the left end of the step, one can hence obtain the solution of linearized SDEs that transversally intersects the true solution at the right end of the step. Here, we exploit the fact that analytical solutions of linear time-invariant SDEs under additive noises are known in closed form. Hence, the quality of solutions via the weak form of stochastic locally transversal linearization (WSLTL) crucially depends on how we conditionally linearize the nonlinear drift field and convert the multiplicative diffusion field into a conditionally additive one. Instead of accounting for drift nonlinearity in a conditionally derivable fundamental solution matrix (FSM) of the linearized system (Roy 2001, 2004), a computationally faster alternative is to treat the nonlinear drift and multiplicative diffusion fields, respectively, as a conditional force vector and a conditionally additive diffusion matrix in the linearized system (Roy *et al*. submitted). We refer to the weak form of the latter procedure as WSLTL-1. Following a somewhat similar logic, here we expand the nonlinear drift and multiplicative diffusion fields conditionally by backward Euler and Newmark schemes. In the process, the replaced nonlinear drift and multiplicative diffusion terms in the linearized system are of higher orders vis-à-vis those in WSLTL-1. The WSLTL schemes via backward Euler and Newmark expansions will be respectively referred to as WSLTL-2 and WSLTL-3 (with the associated strong forms denoted as SLTL-2 and SLTL-3). We also propose simple modifications that allow explicitizations of these implicit schemes. An exposition of the weak replacements of the resulting multiple stochastic integrals (MSIs) is then provided. Local error estimates for the displacement and velocity vectors suggest that the local and global error orders in the new implementations are higher than those in WSLTL-1 and, indeed, the most available schemes. The fact that these higher formal orders also translate into remarkably improved numerical accuracy is demonstrated through implementations of (W)SLTL-2 and (W)SLTL-3 in the context of a few low-dimensional, nonlinear oscillators of general engineering interest.

## 2. Methodology

In the context of engineering systems, a natural starting point would be to consider the following *n* degrees-of-freedom nonlinear oscillator:(2.1)where are the *n*-dimensional state vectors; *C* and *K* are the *n*×*n* state-independent damping and stiffness matrices, respectively; and is the nonlinear part of the (velocity) drift vector. * F*(

*t*) denotes the externally applied deterministic force vector. The diffusion vectors are coefficients of

*r*independently evolving Wiener processes {

*W*

_{r}(

*t*)|

*r*∈

*Z*

^{+},

*r*≤

*q*}, with

*W*

_{r}(0)=0, and

*[*

**E***W*

_{r}(

*t*)

*W*

_{r}(

*s*)]=

*t*∧

*s*. Here,

*denotes the expectation operator. Without a loss of generality, we can also assume that may be decomposed as , where*

**E***σ*

_{r}(

*t*) and are respectively the additive (possibly time-dependent) and multiplicative (state-dependent) diffusion components. Since the Wiener process

*W*

_{r}(

*t*) is continuous and nowhere differentiable in time, (overdot denoting time differentiation) is written only in a formal sense. Thus, we may recast equation (2.1) as the following 2

*n*-dimensional and incremental state space equations:(2.2)where

*A*(

*X*

_{1},

*X*

_{2},

*t*)=−

*KX*

_{1}−

*CX*

_{2}+

*F*(

*t*)−

*A*

_{nl}(

*X*

_{1},

*X*

_{2},

*t*) is the velocity drift vector field written as a function of displacement (

*X*

_{1}≔

**) and velocity () components of the 2**

*X**n*-dimensional response vector . Note that the linear part of the velocity drift field takes the form

*A*

_{l}(

*X*

_{1},

*X*

_{2},

*t*)=−

*CX*

_{2}−

*KX*

_{1}+

*F*(

*t*). Here, we write the displacement and velocity vectors component-wise as and , respectively. Likewise, the vectors

*A*

_{l},

*A*

_{nl},

*σ*

_{r}and

*λ*

_{r}may be written in terms of their scalar components as , , and , respectively. In order to define the solution

*X*(

*t*) properly, we need to use the stochastic integral for a class of right continuous functions

*f*: [

*a*,

*b*]×

*Ω*→

*. We note that this integral is different from the Riemann–Stieltjes integral in the sense that the solution differs depending upon the choice of the point (say left, right or centre) over each time discretization interval where approximating sums of the integrand are taken. The most widely used choices are those due to Ito and Stratonovich, corresponding to the point being located respectively at the left and middle of the interval. While we here interpret the integrals in the Ito sense, it is a straightforward exercise to transform between Ito and Stratonovich integrals (Kloeden & Platen 1999). Now, to ensure that (boundedness of the solution vector*

**R***X*), we require the drift vector, , and the diffusion vector, , to be measurable functions satisfying Lipschitz conditions(2.3)for some positive (real) constant

*Q*and . Here, |.| denotes the Euclidean norm. Moreover, let the initial conditions be mean-square bounded, i.e. and independent of the filtration generated by Wiener processes

*W*

_{r}(

*t*). Thus, the sample continuity of the stochastic flow , for any event

*ω*∈

*Ω*,1 is assured. We order the interval [0,

*T*] of interest as 0=

*t*

_{0}<

*t*

_{1}⋯<

*t*

_{i}⋯<

*t*

_{f}=

*T*and

*h*

_{i}=

*t*

_{i+1}−

*t*

_{i}, where

*i*∈

*N*. Towards a simplified treatment without compromising the generality, a uniform time-step size

*h*

_{i}=

*h*∀

*i*is presently assumed.

The key idea behind (W)SLTL is to replace the nonlinear set of SDEs with an equivalent linearized set, such that the tangent fields to the nonlinear and linearized vector fields remain transversal almost everywhere in the extended state (phase) space and the linearized solution manifold transversally intersects the true solution manifold at the grid points on the time axis. Such transversality arguments enable the construction of the linearized system without using Taylor-like expansions at any stage. Since the integration scheme will be recursively applied over successive time-steps, we consider the incremental form of SDEs (equation (2.2)) over the *i*th discretized time-interval *T*_{i}≔(*t*_{i}, *t*_{i+1}], wherein solution to the *i*th linear system should ‘closely’ represent the nonlinear response over the same interval and must as well transversally intersect the true solution at *t*=*t*_{i+1} while satisfying the initial conditions at *t*=*t*_{i}. Let the initial condition vector over *T*_{i} be , i.e. we assume that is known. Recall that, by way of imposing the transversal intersection, it is necessary to enforce the constraint identity , where denote the transversally linearized WSLTL solution. We are primarily interested in such schemes with higher accuracy and a weak stochastic solution of order ‘*k*’ requires that the following inequality holds for any (sufficiently smooth) deterministic (scalar or vector) function :(2.4)For further work, without a loss of generality, we assume *Ψ* to be a scalar function. In order to obtain the linearized vector field for using the older version of SLTL, we replace the nonlinear part of the drift field and the diffusion field corresponding to the multiplicative noise by using the still-unknown state vector *X*_{i+1} at *t*=*t*_{i+1}. Thus, the WSLTL-based SDEs take the (implicit) form (Roy & Dash 2005*a*)(2.5)This version of weak SLTL will here be referred to as WSLTL-1. Here, conditional replacements of nonlinear drift and multiplicative diffusion fields satisfy the minimal requirement of the vector fields of equations (2.2) and (2.5), remaining instantaneously identical at *t*_{i+1}. No attention is, however, paid to how close the linearized solution is to the true solution. In order to achieve higher orders, we propose to replace *A*_{nl}(*X*_{1}, *X*_{2}, *t*) and *λ*_{r}(*X*_{1}, *X*_{2}, *t*) via implicit substitutions based on backward Euler and Newmark expansions. Towards implicit Euler substitutions of *A*_{nl}(*X*_{1}, *X*_{2}, *t*) and *λ*_{r}(*X*_{1}, *X*_{2}, *t*), the state variables *X*_{1}(*t*) and *X*_{2}(*t*) (appearing as arguments of *A*_{nl} and *λ*_{r}) are backward expanded to *O*(|*t*−*t*_{i+1}|) for *t*∈*T*_{i} as(2.6)

(2.7)In the present weak form of SLTL, we replace the Gaussian stochastic integrals with a set of random variables having appropriately chosen discrete probability distributions that are generated with significantly less computational effort. Hence, the functions *A*_{nl} and *λ*_{r} replaced would be referred to as and , respectively. Note that the required identities and are satisfied. For example, if , then we haveSuch a substitution will result in the following linearized SDEs:(2.8)The weak linearization method based on the above SDEs will henceforth be referred to as WSLTL-2. In order to impose the transversality argument at the right end of sub-interval *T*_{i}, the directional (Gateaux) derivatives of linearized and nonlinear vector fields (in equations (2.2) and (2.8)) should not be identical. This in turn implies that the following inequalities of Frechet derivatives (Jacobians) or must hold for at least one *r*∈[1, *q*] at *t*=*t*_{i+1} (see Roy (2001) for more details on transversality of vector fields).

Note that, as in (W)SLTL-1, (W)SLTL-2 also does not need computation of derivatives of the vector field. In order to achieve a still higher order accuracy without a need to compute derivatives, the functions *A*_{nl}(*X*_{1}, *X*_{2}, *t*) and *λ*_{r}(*X*_{1}, *X*_{2}, *t*) may be conditionally replaced via backward Newmark expansions (Roy & Dash 2002) of *X*_{1} and *X*_{2} over *t*∈*T*_{i}. Note that the displacement and velocity components are expanded to different orders in the Newmark expansion. Thus, while we use the backward Euler expansion (2.7) for *X*_{2} appearing in *A*_{nl}(*X*_{1}, *X*_{2}, *t*) and *λ*_{r}(*X*_{1}, *X*_{2}, *t*), the displacement vector *X*_{1} is expanded to a higher order. The backward expansions for *X*_{1}(*t*) and *X*_{2}(*t*), used in the present form of WSLTL (henceforth called WSLTL-3), are given by(2.9)(2.10)where *α* and *β* are the Newmark implicitness parameters. Numerical experiments have revealed that the (W)SLTL solution is reasonably insensitive to the choice of these parameters and here we use *α*=*β*=0.5 consistently. The conditionally linearized, backward Newmark replacements of *A*_{nl}(*X*_{1}, *X*_{2}, *t*) and *λ*_{r}(*X*_{1}, *X*_{2}, *t*) will be respectively denoted as and .

Upon such replacement of *A*_{nl}(*X*_{1}, *X*_{2}, *t*) and *λ*_{r}(*X*_{1}, *X*_{2}, *t*) via backward expansions, the closed-form solution of the linearized SDEs corresponding to (W)SLTL-2 or (W)SLTL-3 over *T*_{i} may be conditionally written as(2.11)where is the 2*n*-dimensional (suitably pre-augmented with the *n*-dimensional zero vector {0}) force vector and is the similarly augmented 2*n*-dimensional additive diffusion vector. Similarly, the vectors and are suitably augmented (with the zero vector) forms of and , respectively, with the superscript *B* standing for *N* or *E* as appropriate. Moreover, the FSM corresponding to the given SDEs with only the linear drift field is defined as , wherewhere [0]_{n×n} and [*I*]_{n×n} are respectively zero and identity matrices of size *n*×*n*. Finally, the discretized solution vector is found by enforcing the following constraint using equation (2.11) at *t*_{i+1}:(2.12)and then solving the resulting nonlinear algebraic equations for via the Newton–Raphson or fixed-point iteration technique (with the latter being the preferred choice as it involves no derivative calculations). Observe that equation (2.12) is a 2*n*-dimensional system of nonlinear algebraic equations in , provided that is replaced by the r.h.s. of equation (2.11) at *t*=*t*_{i+1}. For WSLTL formulations, we are presently looking for weak replacements of stochastic integrals in equations (2.9) and (2.10), which have the generic form , where is a deterministic function of *s* (the time variable). On expanding about *t*=*t*_{i} via Taylor series, the above integral may be written as(2.13)As we are interested in evaluating the stochastic integral up to a certain order of accuracy (table 1), inclusion of the first few (here not more than four) terms in the Taylor expansion of should suffice. Equation (2.13) may thus be approximated as(2.14)Before going into the details of weak replacements of stochastic integrals, it would be in order to briefly describe the notations presently employed for weak integrals. First, we use and , *Χ*≥1. For *Χ*=1, can be equivalently written as . Similarly, we have . In order to generate realizations of the Wiener component process *W*_{r}(*r*∈[1, *q*]), independent Gaussian random variables are required over each *T*_{i} to evaluate the integrals *ξ*_{r}, , and for SLTL-3. However, we would only need the integrals *ξ*_{r}, and for SLTL-2. For conciseness, we derive the weak distributions only for WSLTL-3, with the understanding that these arguments are readily adaptable for WSLTL-2. Apart from *ξ*_{r} and , there would be additional random variables owing to backward Euler and Newmark substitutions of the nonlinear drift and multiplicative diffusion terms. For instance, using backward Euler expansion, integrals of the form , *Χ*≥1 arise. Now a weak replacement of the r.h.s. of equation (2.14) is written as (here *ϑ*_{r} and denote weak replacements of *ξ*_{r} and , respectively)(2.15)The above equality is only in terms of the associated probability measures. Similarly, let be the weak replacement for the additional random variable (stochastic integral generated owing to the backward expansion). We may readily check that are uncorrelated with , where *r*, *m*∈[1, *q*] and (Kloeden & Platen 1999). Note that the weak variables *ϑ*_{r}*h*^{0.5} and are of respective orders *O*(*h*^{0.5}) and *O*(*h*^{i−0.5}). A formal analysis of the orders of such weak approximations is provided in §3. However, here we note that (*δ*_{rm} is the Kronecker delta) and that the expected value of *ξ*_{r}, or with is at least of order *O*(*h*^{4.0})—which is equal to or higher than the desired order of (local) error for computing the expected values of functions of response states using WSLTL-3. Thus, it suffices to model deterministically as (with probability 1). Accordingly, for weak solutions, we require to characterize only the weak variables *ϑ*_{r} and , presently with discrete and symmetric distributions about the zero mean, for each *r*∈[1, *q*], such that all the associated odd moments (up to *O*(*h*^{3.5})) are zero. In particular, we have(2.16a)(2.16b)(2.16c)(2.16d)(2.16e)(2.16f)(2.16g)We now need to specify a symmetric joint distribution of *ϑ*_{r} and *ϑ*_{1r} about the origin (*ϑ*_{r}, *ϑ*_{1r})=(0, 0). Towards this, we need to determine the even-order moments , , , and , *r*, *m*∈[1, *q*]. It is observed that all other even-order moments among different elements of {*ϑ*_{r}*h*^{0.5}| *r*=1, …, *q*}, and have orders higher than *O*(*h*^{4}), and thus they need not be accounted in the weak approximation. For modelling , we just need a two-point distribution (symmetric about the origin), with variance . At this stage, the following relations may be readily derived (see §3 for more details):(2.17a)

(2.17b)

(2.17c)

(2.17d)(2.17e)(2.17f)As modelling of Gaussian random variables requires the usage of transcendental functions (sine, cosine and logarithm), their weak replacement through discrete distributions speeds up the computation manifold. Presently, we derive the joint probability distribution of the random variables required in WSLTL-3. Here, we need to determine the joint probability distribution of *ϑ*_{r}*h*^{0.5}, and satisfying the following expectations:(2.18)Additionally, discrete distributions of and must satisfy(2.19)We need to find out the following quantities: *q*_{0}, *q*_{1}, *q*_{2}, *q*_{3}, *q*_{4}, *k*, *l*, *m*, *r*_{1}, *r*_{2}, *r*_{3}, *r*_{4}, *r*_{5}, and *r*_{6}. Any symbolic manipulator may be used for this purpose and we have presently used Mathematica. The following resulting solutions are:(2.20)Hence, the joint distribution of weak random variables appearing in WSLTL-3 is given by(2.21)Similarly, we can also obtain the following joint distribution for the weak random variables corresponding to WSLTL-2:(2.22)

### (a) Explicit forms of the present schemes

In the above implicit forms, we conditionally replace and in equation (2.2) with their backward Euler equations (2.6) and (2.7) and Newmark expansions (2.9) and (2.10), such that the linearized solution over the interval *T*_{i} would closely represent the nonlinear response and transversally intersect the true solution at *t*=*t*_{i+1}. Now, rather than replacing the nonlinear drift and multiplicative diffusion fields by backward expansions, if we replace them with forward expansions based at *t*=*t*_{i}, linearized solutions may be written explicitly as functions of the known initial condition and they would still be of the same order as achievable through the implicit forms. Such an explicitization could be important, as explicit methods are computationally more efficient. However, as forward expansions of and would be based on , the linearized vector field would not be instantaneously identical with the nonlinear one at *t*=*t*_{i+1} and the transversality argument (see §2, following equation (2.3)) would no longer hold. We now write the forward expansions for *X*_{1}(*t*) and *X*_{2}(*t*) appearing in the arguments of and for explicit Euler substitutions in *t*∈*T*_{i} as(2.23)(2.24)Hence, the arguments *X*_{1}(*t*) and *X*_{2}(*t*) substituted in functions and would be referred to as and , respectively. Similarly, one can use explicit Newmark expansions of *X*_{1}(*t*) and *X*_{2}(*t*) and thereby obtain the explicitly computable linearized functions and . However, while such explicit linearizations would considerably ease the computational overhead and would invariably be the preferred choice in real-time applications, an implicit method may still work better for higher step sizes. Note that modelling stochastic integrals and their weak replacements in the explicit forms follow similar steps as in the implicit forms.

## 3. Local error analysis

In this section, we determine the weak (local) orders of WSLTL-2 and WSLTL-3 over the time-interval (*t*_{i}, *t*_{i+1}]. Error estimates for the weak forms are somewhat similar to those derived for the higher order strong forms of SLTL (Saha & Roy 2007). The orders of local errors are accordingly governed by those of the leading-order terms that do not agree with each other in stochastic Taylor expansions of conditionally linearized and true solutions over *T*_{i} using the corresponding vector fields. In order to prove the weak order of convergence, some additional work is however necessary to account for the weak replacements of the Gaussian stochastic integrals by those with a discrete distribution.

### (a) The Ito–Taylor expansion

Let *Ψ*(*X*_{1}, *X*_{2}, *t*) (assumed to be a scalar function without a loss of generality) be continuous and sufficiently smooth in its arguments, with *X*_{1}, *X*_{2} being governed by the SDE (2.2). Moreover, let *Ψ*(*X*_{1}, *X*_{2}, *t*) be adapted to the filtration generated by *W*_{r}(*t*). Then, *Ψ*(*X*_{1}, *X*_{2}, *t*) is also an Ito process whose SDE may be written as(3.1)The operators *Λ*_{r} and *L* are given by(3.2)(3.3)The basis for the Ito–Taylor expansion is to expand *Ψ*(*X*_{1}, *X*_{2}, *t*) as(3.4)Higher order expansions of this basic expansion are readily achievable by repeatedly applying Ito's formula to functions appearing within the integrals (such as *Λ*_{r}*Ψ* and *LΨ* in the r.h.s. of equation (3.4)). We refer to Kloeden & Platen (1999) for further details.

A MSI has the typical form(3.5)Here, *j*_{1}, *j*_{2}, …, *j*_{k} are integers taking values in the set {0, 1, …, *q*} and is referred to as the *k*th (Ito) MSI. Moreover, d*W*_{0}(*s*) is taken to indicate d*s*.

A function belongs to the class *C*_{ρ} if(3.6)We assume that *C*_{ρ} contains all the integrands in the Ito–Taylor expansion (including the remainder). Equivalently, we assume that the drift coefficient *a*, the diffusion coefficient *b*_{r} (2.3) and their partial derivatives of sufficiently high order are in *C*_{ρ}.

We use the notations , and . Moreover, the triple and double integrals over [*t*_{i}, *t*_{i+1}] are written as and , respectively.

### (b) Error analysis for weak form of stochastic locally transversal linearization-2

For brevity, we provide the error analysis for WSLTL-2. The analysis for WSLTL-3 is analogous and we thus state only the final error orders. We need to compare the Ito–Taylor expansions of linearized and true solutions (i.e. the corresponding displacement and velocity components) over *T*_{i} in terms of their associated vector fields based on the same initial conditions at *t*_{i}. The terms that do not match with each other in these expansions presently constitute the remainder terms, whose leading-order terms determine the local error orders. For conciseness, we derive such remainder terms for displacement only and just state those for velocity.

### (c) Remainder terms for displacement and velocity

Let the displacement component , corresponding to equation (2.2), be expanded as(3.7)In WSLTL-2, we replace with (superscript *E* denoting backward Euler approximations of the argument states), which is given by (in the vector form)(3.8)where *t*∈*T*_{i}, is the Frechet derivative (matrix) of the vector *A*_{nl} at . Moreover, is the remainder vector corresponding to a one-step approximation of and is of order *O*(*t*−*t*_{i+1}). In the specific case of system (2.2) being driven only by additive noises, will be *O*(*t*−*t*_{i+1})^{1.5}. Towards estimating the error due to the backward Euler replacement, we expand in through a backward Ito–Taylor series as(3.9)The remainder is given by(3.10)Now expanding (r.h.s. of equation (3.9)), in a generalized Taylor expansion and using equation (3.8), we get(3.11)where contains higher order derivatives of the Taylor expansion. In general, is of order *O*(*t*−*t*_{i+1})^{1.5} and in case of the system (2.2) being driven by only additive noises, is of order *O*(*t*−*t*_{i+1})^{2}. Hence, the signed error vector due to backward Euler replacement of is expressible as(3.12a)The order of is governed by , as the other terms are of still higher orders. Consequently, while is generally of order *O*(*t*−*t*_{i+1}), the order of accuracy will increase to *O*(*t*−*t*_{i+1})^{1.5} when the system (2.1) is driven by additive noises only.

Now, a similar exercise is needed to account for the replacement of by . As in the case of , remainder term for the multiplicative diffusion case would be of order *O*(*t*−*t*_{i+1}). Moreover, the other remainder terms and would also be *O*(*t*−*t*_{i+1}). Hence, the signed error due to backward Euler replacement of is(3.12b)On replacing and via backward Euler approximations in equation (3.7), we finally get(3.13)Hence, the error component is of the following form due to linearization (over *T*_{i}):(3.14)Likewise, we expand the velocity component of equation (2.2) as(3.15)with and replaced by and , respectively. Recall that and are of order *O*(*t*−*t*_{i+1}) (under combined multiplicative and additive noises) or of order *O*(*t*−*t*_{i+1})^{1.5}(under only additive noises). Equation (3.15) is now written as(3.16)Here, the local error component for due to linearization is of the form(3.17)

### (d) Strong error orders

We state the following bounds on expectations of remainder terms and their products with stochastic integrals associated with the displacement and velocity for SLTL-2.

*Under conditions following equation* (3.6), *Lipschitz continuity* (2.3), *one has the following estimates for* *and* *j*∈[1, *n*], *in SLTL*-2:(3.18)(3.19)(3.20)(3.21)*for some scalar function* *independent of h*.

From the expression of *R*_{1} in equation (3.14), it follows that (as expectations of the terms associated with would be identically zero)(3.22)Since and (3.6) for every *j*∈[1, *n*], one can find an even integer 2*ζ*, such that(3.23)Substituting (3.23) into (3.22) and noting that is bounded by some quantity (Milstein 1995), the first inequality of (3.18) follows. In a similar way, the other inequalities of (3.18)–(3.21) may be shown (Roy 2006). ▪

Now, introduce the following notations for *n*-dimensional exact (or true), strong and weak response increment vectors as:(3.24)The left superscript *S* is used to denote (increments of) the response variables via the strong SLTL. While *Δ*_{1} and *Δ*_{2} are the true response increments, ^{S}*Δ*_{1} and ^{S}*Δ*_{2} are those via the strong SLTL and and are those via the WSLTL. One has the following relations between strong and exact increments:(3.25)The following proposition establishes the ‘closeness’ (in terms of orders of *h*) of the first few moments of various elements of true and strong increment vectors.

*For Lipschitz bounded drift and diffusion vectors* (2.3), *the following inequalities hold*:(3.26)(3.27)(3.28)*for j*_{k}, *j*_{p}∈*Z*^{+} *and some scalar functions* . *Further, one has*(3.29)

From equation (3.25) and given that local initial conditions are being treated as deterministic, it follows that (via proposition 3.1)To prove the inequality (3.26) for *m*=2 (i.e. for the second moments of displacement increments), it is first noted (via equation (3.25)) that(3.30)Now, using equations (3.14) and (3.18), one has , and also it follows from the expression for ^{S}*Δ*_{1} that (via Cauchy–Schwarz inequality). Thus, the r.h.s. of equation (3.30) is ≤*O*(*h*^{2}). Similarly, inequality (3.27) for *m*=3 may also be proved. Inequality equation (3.27) for *m*=1 is also readily proved, given that . Then, for *m*=2, one has(3.31)Referring to equation (3.17), one has , and it is only to be shown that the other two terms on the r.h.s. of equation (3.31) are also at least *O*(*h*^{2}). Towards this, some of the lowest order terms of *R*_{2} (i.e. those of orders *O*(*h*^{1.5}) and *O*(*h*^{2})) are Ito–Taylor expanded and then multiplied with the Ito–Taylor expansions of the (strong) SLTL-based velocity components to yield(3.32)for some . From the first principles of Ito calculus, it may be shown that . Thus, we may write(3.33)We refer to Milstein (1995) or Roy (2006) for the derivation of the above expression. The rest of the inequalities (3.27)–(3.29) may also be shown in a similar manner.

### (e) Weak error orders

Consider the WSLTL-2 scheme (2.8)–(2.10). Using the weak random variables , the WSLTL-2 displacement map is written as (*p*∈[1, *n*])(3.34)

Equivalently, the WLTL-2-based displacement increment is given by(3.35), , and are given by(3.36)Similar expressions for weak velocity components may also be written. Note that the expression for strong displacement increment, , is given by equation (3.35), with replaced by . For oscillators under only additive excitations, equations (3.35) and (3.36) provide an approximation with a local order of accuracy *O*(*h*^{2.5}) for the displacement components (or increments), provided that the weak replacement of random variables is adequately accurate. A similar observation is valid for the velocity components, with the exception that they are locally approximated to order *O*(*h*^{1.5}). Statistical properties of *ϑ*_{r}*h*^{0.5}, , and had earlier been derived in equations (2.16*a*)–(2.18).

The following proposition establishes the local order of accuracy of WSLTL-2.

*Let the Lipschitz continuity requirements as in equation* (2.3) *hold*. *Also suppose that the moment relations in equations* (2.16a)–(2.16g) *and* (2.17a)–(2.17f) *hold*. *Then, one has*(3.37)(3.38)(3.39)*for i*_{j}, *i*_{k}∈*Z*^{+} *and* , , *independent of h. Further, one has*(3.40)(3.41)

Assertions (3.37) and (3.38) follow directly from the WSLTL-based expressions for and (3.35) by noting that and have terms containing *ϑ*_{r}*h*^{0.5}, , and . It is also noted from Ito–Taylor expansions for and (based respectively at *X*_{1,i} and *X*_{2,i}), whose evolutions are governed by the WSLTL-2-based SDEs (2.5), that the least order terms in and are respectively of orders *O*(*h*) and . It is now simple to see (via direct computations, using inequalities (3.18)–(3.21) along with Cauchy–Schwarz inequality) that if the moment relations of equations (2.16*a*)–(2.16*g*) and (2.17*a*)–(2.17*f*) hold, then inequalities (3.39)–(3.41) also hold. Note that, under additive noises only, l.h.s. of inequalities (3.37) and (3.40) would be bounded by and l.h.s. of inequalities (3.38), (3.39) and (3.41) would be bounded by or as applicable. Together with inequalities (3.26)–(3.29) of proposition 3.2, this means that at least first three moments of strong and weak LTL-based increments of displacements alone are of the same error order. Moreover, first four moments of strong and weak LTL-based velocity increments and first three mixed moments of pathwise and weak LTL-based displacement and velocity increments are also of similar error orders.

Based on the preceding discussion, the main result on the error orders of the WSLTL-2 method may now be proposed in the form of the following theorem.

*Let all the conditions as stated in* *propositions* 3.1–3.3 *hold*. *Also suppose that Ψ*_{1}(*X*_{1}, *t*)∈*C*_{ρ} *is a function of displacement alone and Ψ*_{2}(*X*_{1}, *X*_{2}, *t*)∈*C*_{ρ} *is a function of both displacement and velocity vectors. Moreover, assume that all partial derivatives of Ψ*_{1} *with respect to X*_{1} *up to order 4 and all partial derivatives of Ψ*_{2} *with respect to X*_{1} *and X*_{2} *up to order 5 exist and they belong to C*_{ρ} *as well*. *Then, one has, for the WSLTL*-2 *approximations* , *the following inequalities*:(3.42)(3.43)

Note that and , for *m*=3 and *p*+*r*=3. Now, the proof for the inequality (3.42) may proceed on the following lines (Milstein 1995). Expand *Ψ*_{1}(*X*_{1,i+1})=*Ψ*_{1}(*X*_{1,i}+Δ_{1}) via a Taylor expansion in integer powers of the displacement increment vector *Δ*_{1}=*X*_{1,i+1}−*X*_{1,i} based at *X*_{1,i} and leaving out Lagrange's remainder terms of order *h*^{3} and higher. Similarly, expand in powers of the weak displacement increment vector based at *X*_{1,i} and again leaving out Lagrange's remainder terms of order *h*^{3} and higher. Now, inequality (3.42) follows if one takes the norm of the difference of expectations of these expansions followed by the use of propositions 3.2 and 3.3. Inequality (3.43) may also be derived in a similar manner. When the WSLTL-2-based SDEs are subjected to additive noises only (i.e. ), then the l.h.s. of inequalities (3.42) and (3.43) would be bounded by and , respectively. ▪

### (f) Distributions of weakly replaced stochastic integrals

Finally, the following proposition is used to justify the simplified discrete distributions of the weak random variables *ϑ*_{r}*h*^{0.5}, , and .

*The moment relations in equations* (2.16a)–(2.16g) *and* (2.17a)–(2.17f) *are satisfied by a family of discrete or continuous distributions of ϑ*_{r}*h*^{0.5}, , and . *A convenient (symmetric and discrete) distribution is given by equation* (2.18) *that is reproduced in the following for ready reference*:(3.44)

First, it is noted that , which is a weak replacement for the Gaussian random variable , can be directly modelled using the random variable along with a deterministic multiplier. This follows since, except for , all other moments of with the other two random variables and with itself are either zero or of orders higher than *O*(*h*^{3}), for all *r*, *m*∈[1, *q*]. Now consider the joint distribution of and for any *r*. If a symmetric (and discrete) distribution around is chosen, then it ensures the vanishing of all odd-order moments involving . Similarly, as is uncorrelated with , and as we need to satisfy up to second moment of (i.e. ), symmetric distribution about 0 with equal probability masses on either side would suffice. Let this distribution be given by the last of equation (3.44). Then, it is simple to verify that all the required moment equalities in equations (2.18) and (2.19) are identically satisfied. By actual construction, it may be shown that such a discrete distribution is not unique. For instance, the following is also a valid distribution for :(3.45)This completes the proof of the proposition.

The analysis of local error orders for WSLTL-3 constitutes a similar exercise. We briefly report the global error order for different schemes in table 1.

Finally, we note that weak error orders for the explicit linearization schemes (§2) remain the same as those for the implicit WSLTL schemes, provided that distributions of the weak random variables are appropriately derived.

## 4. Numerical illustrations

In order to demonstrate the performance of WSLTL-2 and WSLTL-3 vis-à-vis WSLTL-1, we numerically simulate strong and weak responses of Duffing and Van der Pol oscillators under additive and/or multiplicative noises.

### (a) The Duffing oscillator

Consider the five-parameter hardening Duffing equation in the incremental form(4.1)On linearized replacements of the nonlinear drift and multiplicative diffusion terms (depending on the form of transversal linearization, see §2), the SLTL-based replacement of equation (4.1) is given by(4.2)Recall that, depending on the form of WSLTL, the superscript *B* is replaced by *E* (WSLTL-2) or *N* (WSLTL-3). The linearized solution can be written conditionally over *t*∈*T*_{i} as(4.3)where is the FSM corresponding to the linear part of the drift field,, , and . For clarity, we write the following explicit forms of :(4.4)

(4.5)Similar expansions may be performed to obtain and . We validate some of our results with available ‘exact stationary solutions’ of the associated probability density functions (pdfs) of the reduced Fokker–Planck equation (Wang & Zhang 2000). For the present oscillator under only additive noises, the exact stationary pdf is given by(4.6)where the constant *A* has to be found so as to satisfy the normalization constraint . Presently, we have used the symbolic manipulator Maple for finding the expectation of any deterministic function as . The numerical integration used is an adaptive Newton–Cotes rule with a sufficiently small tolerance (as compared with the error orders of the methods). For non-removable endpoint singularities, an adaptive double-exponential quadrature method has been used.

In figure 1, we provide the history plots of * E*[

*X*

^{2}] and only under a low additive noise (small

*ϵ*

_{4}) for different time-step sizes using WSLTL-1, WSLTL-2 and WSLTL-3 (with a consistent ensemble size of 5000). We also compare the simulation results with the exact stationary solution for large

*t*. For

*h*=0.1, while solutions through WSLTL-1 overflow, we obtain numerically stable (even though somewhat inaccurate) solution histories using WSLTL-2 and WSLTL-3 (figure 1

*a*). However, as shown in figure 1

*b*, a reduction in the step size to

*h*=0.01 enables all the three different WSLTL strategies to remain close to the correct stationary solution for sufficiently large

*t*. Figure 2 shows that, under high additive noise (large

*ϵ*

_{4}), results via WSLTL-2 and WSLTL-3 remain close to their strong counterparts. However, while the results via WSLTL-3 (and hence strong SLTL-3) remain close to the exact stationary solution for large

*t*(figure 2

*b*), WSLTL-2 (and strong SLTL-2) yield results with conspicuous inaccuracy (figure 2

*a*). This is a substantive demonstration of the importance of higher order schemes for stochastic dynamical systems. Attempts to use WSLTL-1 for this case (with

*h*=0.01) lead to solution overflows. In table 2, we tabulate the consumed CPU times for the three variants of the SLTL method. The results emphasize the computational efficiency of the WSLTL schemes over the strong SLTL schemes. All the programs have presently been run on Compaq Alpha ES40 systems under a similar load-sharing facility.

We now consider the Duffing oscillator under combined sinusoidal force (deterministic) and additive noise such that chaotic behaviour is exhibited by the oscillator for certain parameter choices. In figure 3, we plot the phase portraits using the three variants of WSLTL for step sizes of *h*=0.001 and *h*=0.01. While performance of WSLTL-3 remains nearly invariant, the quality of solution via WSLTL-2 deteriorates and via WSLTL-1 again overflows for *h*=0.01. In figure 4, we show time history plots of second moments under combined deterministic excitation and high multiplicative noise. Once again, we find that present versions of WSLTL have remarkably improved performance over WSLTL-1. The numerically superior and computationally stable performance of WSLTL-3 over its two other counterparts is once more evidenced. Though not explicitly shown in these figures, time histories via WSLTL-1 blows up while those through WSLTL-2 yields computationally stable yet numerically inaccurate solutions for larger step sizes of *h*=0.01. On the other hand, the WSLTL-3 works fine with *h*=0.01.

### (b) The Duffing–Holmes oscillator

We now consider a double-well Duffing–Holmes oscillator subjected to additive noise only. The governing SDEs are(4.7)If the deterministic force is absent and the oscillator is driven by only an additive noise with low intensity, then, for large *t*, we may expect to have a stationary bimodal pdf with the modes around the couple of (stable) fixed points (corresponding to the unforced, damped oscillator) symmetrically placed about the saddle at (*X*_{1}, *X*_{2})=(0, 0) on the *X*_{1}-axis. Figure 5 shows snapshot plots of pdfs obtained using WSLTL-3 in the transient regime and as stationarity is approached. The nature of the pdfs obtained in figure 5 is similar to that reported by Kunert & Pfeiffer (1991). Interestingly, the small hump at the saddle (i.e. at (*X*_{1}, *X*_{2})=(0, 0)) indicates an onset of diffusion across the separatrix. Note that simulations using 50 000 samples are carried out for obtaining the pdfs. Though better methods are available (Docenko & Berzins 2003), the pdfs are constructed by the standard procedure of binning the space and approximating the pdf by the ratio of the number of events falling inside each bin over the total and normalized to the bin volume. A cubic spline smoothing is then adopted before reporting the final pdfs.

### (c) The Van der Pol oscillator

Unlike the Duffing oscillator, the Van der Pol oscillator exhibits nonlinear damping. The governing equation, under additive noise, is of the form(4.8)Without getting into a more detailed reporting of such results, we note that the WSLTL-3 and WSLTL-2 continue to show numerically superior performance over WSLT-1 for this oscillator too. In particular, WSLTL-1 consistently exhibits oscillations with considerable larger amplitudes than WSLTL-2 and WSLTL-3, especially for larger step sizes. Moreover, with increase in the nonlinearity parameter *μ*, spurious oscillations exhibited by the lower order methods (i.e. WSLTL-1 and WSLTL-2) are sharply reduced in simulations via WSLTL-3. Figure 6 shows plots of stationary pdfs obtained using WSLTL-3 under different values of *μ*. The nature of pdfs obtained through the WSLTL-3 maps is similar to those reported in the literature (Naess & Haestad 1994).

## 5. Conclusions

This work deals with the issue of deriving accurate and efficient numerical integrators for stochastically driven nonlinear oscillators. While numerical accuracy is important given a paucity of available higher order schemes for SDEs, computational efficiency is crucial in Monte Carlo simulations of large engineering systems as well as in real-time estimation problems. We have proposed two derivative-free modifications (referred to as SLTL-2 and SLTL-3) of the LTL method and thus obtained solutions with significantly enhanced numerical accuracy. This is possible through linearized replacements of the nonlinear drift and multiplicative diffusion terms via backward Euler or Newmark expansions over a time-step. Such backward expansions enable the linearized fields to remain identical with the original ones at the forward time point (of a given time-step) and also ensure a higher closeness between the two vector fields. This forms the basis of arguments that indicate higher formal orders achievable through the present schemes vis-à-vis SLTL-1 and, indeed, most existing schemes in the literature (Burrage *et al*. 2004; Roy 2006; Komori 2007). Apart from the derivative-free nature of this class of linearizations, the other source of simplicity is in the usage of not-too-many stochastic integrals, especially for the weak formulations. These stochastic integrals are originally Gaussian in nature. However, as emphasized in this work, their weak modelling through discrete probability distributions speeds up the computation 2.5–3 times. Unlike strong schemes, an added simplification in such weak formulations is that a few of the higher order stochastic integrals become identically zero. We must emphasize that the strength of the present linearization is in its simplicity of implementation irrespective of the dimensionality of the oscillator.

## Acknowledgments

The authors acknowledge the financial support from the Department of Atomic Energy, Government of India, through contract no. DAEO/MCV/DBR/112.

## Footnotes

↵† Graduate student.

↵

*Ω*is the sample space, defined on the probability space of [*Ω*, , ] of which the pair (*Ω*, ) are measurable space with a total measure of 1, i.e. [*Ω*]=1 and :→[0, 1] is a probability measure transforming under Kolmogorov axioms.- Received January 31, 2007.
- Accepted April 5, 2007.

- © 2007 The Royal Society