First Passage Time of Nonlinear Diffusion Processes with Singular Boundary Behavior

TL;DR

This paper develops new theorems for calculating the moments of first passage times in nonlinear stochastic processes with singular boundaries, validated through existing results and applied to ship roll dynamics.

Contribution

It introduces analytical theorems for moments of first passage times in nonlinear processes with singular boundaries, applicable to physical systems like ship dynamics.

Findings

Validated theorems using existing analytical results.

Provided fast computational methods for moments of first passage times.

Applied results to model dangerous ship roll behavior.

Abstract

New theorems for the moments of the first passage time of one dimensional nonlinear stochastic processes with an entrance boundary are formulated. This important class of one dimensional stochastic processes results among others from approximations of the energy or amplitude of second order nonlinear stochastic differential equations. Since the diffusion of a stochastic process vanishes at an entrance boundary, the entrance boundary is called a singular point of the stochastic process. The theorems for the moments of the first passage times are validated based on existing analytical results. In addition, the first passage times of a nonlinear stochastic differential equation, which is important for the determination of dangerous ship roll dynamics, are calculated. The proposed analytical expressions for the moments of the first passage times can be calculated very fast using standard quadrature formulas.

Tables

Table 1. Table 1: Parameter values used for computations

${\alpha }_1=3.187$	${\alpha }_3=4.164$	${\beta }_1=0.655$	${\beta }_2=0.921$	${\beta }_3=0$
${\nu }_1=0.018$	${\nu }_2=1.783$	$H_c=0.529$	$b_c={40}^{\circ }$	$\epsilon =0.1$

Equations

$$\mathrm{d}x(t)=m(x(t),t)\mathrm{d}t+\sigma (x(t),t)\mathrm{d}W,$$

$$\mathcal{L}p(x)=m(x)\frac{\partial }{\partial x}p(x)+\frac{1}{2}{\sigma }^2(x)\frac{{\partial }^2}{\partial x^2}p(x)$$

$${\mathcal{L}}^{\ast }p(x)=-\frac{\partial }{\partial x}\left(m(x)p(H)\right)+\frac{1}{2}\frac{{\partial }^2}{\partial x^2}\left({\sigma }^2(x)p(x)\right).$$

$$S(x)={\int }^{x}s(x)dx,$$

$$s(y)=\exp \left(-2{\int }^y\frac{m(\theta )}{{\sigma }^2(\theta )}d\theta \right)$$

$$\mu (y)=1/({\sigma }^2(y)s(y))$$

$$dM(y)=\mu (y)dy,$$

$$\mathcal{L}p(x)=\frac{1}{2}\frac{\partial }{\partial M}\left(\frac{\partial p(x)}{\partial S}\right),$$

$${\mathcal{L}}^{\ast }{p}_{st}(x)=0,$$

$$p_{st}(x)=\mu (x)[c_{1}S(x)+c_2].$$

$${\Sigma }_l(x_l)=\infty \mathrm{and}{N}_l(x_l)<\infty ,$$

$${\Sigma }_l(x_l)<\infty ,N_l(x_l)<\infty \mathrm{and}M(x_l)=0,$$

$${\Sigma }_l(x_l)<\infty \mathrm{und}M(x_l,x]=\infty .$$

$$N_l(x_l)={\int }_{x_l}^{x_0}\left\{{\int }_z^{x_0}s(y)dy\right\}\mu (z)dz={\int }_{x_l}^{x_0}S[z,x_0]\mu (z)\mathrm{d}z,$$

$${\Sigma }_l(x_l)={\int }_{x_l}^{x_0}\left\{{\int }_{x_l}^{z}s(y)dy\right\}\mu (z)dz={\int }_{x_l}^{x_0}S(x_l,z]\mu (z)\mathrm{d}z,$$

$$S(x_l,z]:=\underset{a\downarrow x_l}{\lim }S[a,z],$$

$$S[a,b]:={\int }_a^{b}s(y)\mathrm{d}y,a,b\in \text{mathds}R.$$

$$p_{st}(x)=\frac{c_2}{{\sigma }^2(x)}\exp \left(2{\int }_{x_e}^x\frac{m(x)}{{\sigma }^2(x)}\right).$$

$$\mathcal{L}{M}_n(x_0)=-n{M}_{n-1}(x_0),M_n(x_c)=0,$$

$$M_1(x_0)=2\left[1-\frac{S[x_0,x_c]}{S[△,x_c]}\right]{\int }_{x_0}^{x_c}S[z,x_c]\mu (z)\mathrm{d}z+2\frac{S[x_0,x_c]}{S[△,x_c]}{\int }_{△}^{x_0}S[△,z]\mu (z)\mathrm{d}z.$$

$$M_1(x_0)=2{\int }_{x_0}^{x_c}\left[{\int }_z^{x_c}s(y)\mathrm{d}y\right]\mu (z)dz+2{\int }_{x_0}^{x_c}s(y)\mathrm{d}y{\int }_{x_l}^{x_0}\mu (z)dz$$

$$\underset{△\downarrow x_l}{\lim }{\int }_{△}^z\mu (\eta )\mathrm{d}\eta <\infty ,\forall x_l<z\le x_c$$

$${\Sigma }_l(x_l)+N_l(x_l)=\underset{△\downarrow x_l}{\lim }\{S[△,z]{\int }_{△}^z\mu (\eta )\mathrm{d}\eta \}$$

$$\underset{△\downarrow x_l}{\lim }S[△,z]=S(x_l,z]=\infty ,\forall x_l<z\le x_c.$$

$$\underset{△\downarrow x_l}{\lim }2\left[1-\frac{S[x_0,x_c]}{S[△,x_c]}\right]{\int }_{x_0}^{x_c}S[z,x_c]\mu (z)\mathrm{d}z=2{\int }_{x_0}^{x_c}S[z,x_c]\mu (z)dz.$$

$${\Sigma }_l(x_l)=\underset{△\downarrow x_l}{\lim }{\int }_{△}^{x_0}S[△,z]\mu (z)\mathrm{d}z=\infty .$$

$$Z:=$$

$$Z$$

$$=2S[x_0,x_c]{\int }_{x_l}^{x_0}\underset{△\downarrow x_l}{\lim }\frac{S[△,z]}{S[△,z]+S[z,x_c]}\mu (z)\mathrm{d}z$$

$$=2S[x_0,x_c]{\int }_{x_l}^{x_0}\underset{△\downarrow x_l}{\lim }\frac{1}{1+\frac{S[z,x_c]}{S[△,z]}}\mu (z)\mathrm{d}z.$$

$$Z=2S[x_0,x_c]{\int }_{x_l}^{x_0}\mu (z)\mathrm{d}z,$$

Full text

Keywords: first passage time; random excitation; stochastic dynamics, Duffing oscillator

1 Introduction

A very important problem in the field of diffusion processes and random dynamical systems deals with the determination of times of the considered nonlinear diffusion process to reach certain boundaries for the first time, starting at a prescribed initial condition.

Using approximation techniques such as stochastic averaging, as described for example in Sri Namachchivaya (1991); Roberts (1978); Dostal et al. (2012), many nonlinear problems of higher order can be approximated by one dimensional Itô stochastic differential equations for the total energy or the response amplitude. Such one dimensional stochastic processes have singular boundaries where the diffusion becomes zero.

This makes it necessary to analyze the limiting behavior of these diffusion processes at such boundaries, since standard formulas or numerical calculations can not be applied appropriately without considering these singularities. One dimensional diffusion processes of this kind arise from approximations of the energy or amplitude of important second order nonlinear problems, for example dynamical systems which undergo a co-dimension two bifurcation Sri Namachchivaya (1991), oscillators with nonlinear damping Roberts (1978), pendulum or ship roll dynamics Roberts & Vasta (2000); Dostal et al. (2012); Dostal & Kreuzer (2016); Dostal et al. (2018).

In many applications it is important to make sure that the response process of a system does not leave a safe domain during its operation time. Then the first time the response process leaves the safe domain has to be determined. Being able to establish formulas for such first passage times, leads to a major impact on the determination of the system reliability and on the expected duration before system failures occur.

After introducing the necessary theory of diffusion processes, we state our new theorems for the moments of the first passage time in section 3, followed by three example applications of the presented theory in sections 5 to 7.

In these sections the first passage time of a Linear system with external excitation is analyzed first, in order to illustrate the theory. The obtained results are validated using previously obtained results by Ariaratnam & Tam (1976). Then results of the presented theory for a forced and damped Mathieu oscillator are analyzed and validated. The Mathieu oscillator is widely used for the modeling of parametrically excited physical systems and has been the subject of many studies, for example Ariaratnam & Tam (1976) Vanvinckenroye & Denoel (2017). Finally, results are shown for a Duffing oscillator forced by external and parametric excitation. This nonlinear system has many applications. It is for example important for the dynamics and capsizing analysis of ships Dostal et al. (2012) and for the analysis of energy harvesting Yurchenko et al. (2013); Dostal et al. (2018).

2 One dimensional diffusion process

For many problems it is necessary to consider a time homogeneous regular diffusion $\{x(t),t>0\}$ on the interval $D=(x_{l};x_{c}]$, which satisfies the Itô stochastic differential equation

[TABLE]

where $m:\mathds{R}\times\mathds{R}\rightarrow\mathds{R},$ and $\sigma:\mathds{R}\times\mathds{R}\rightarrow\mathds{R}$ are measurable functions, see e.g. Oksendal (1992).

As shown in the previous section, such processes result for example from stochastic averaging of general nonlinear oscillators.

The differential generator $\mathcal{L}$ of the diffusion $x$ defined by the Itô equation (1) is given by

[TABLE]

and the associated adjoint operator is

[TABLE]

We introduce the scale measure

[TABLE]

where

[TABLE]

is the scale density. After defining the speed density

[TABLE]

and the speed measure $M$ by

[TABLE]

we can transform the operator $\mathcal{L}$ to

[TABLE]

such that the drift is identically zero. A discussion on the meaning of scale density $s$ and speed density $\mu$ can be found in Karlin & Taylor (1981).

The stationary probability density function $p_{st}$ associated with (1), is the solution of the Fokker-Planck equation

[TABLE]

where $\mathcal{L}^{*}$ is the adjoint of the operator $\mathcal{L}$. The solution of (8) can be given in terms of

[TABLE]

The coefficients $c_{1}$ and $c_{2}$ are determined by the boundary and normality conditions.

Following (Karlin & Taylor, 1981) the (left) boundaries can be classified as

• •

Entrance, if

[TABLE]

• •

reflecting, if

[TABLE]

• •

exit, if

[TABLE]

Here, $\Sigma_{l}(x_{l})$ is the time to reach the left boundary $x_{l}$ starting at $x_{0}\in[x_{l},x_{c}]$, whereas $N_{l}(x_{l})$ is the time to reach $x_{0}\in(x_{l},x_{c}]$ starting at $x_{l}$. These measures are given by

[TABLE]

and

[TABLE]

where

[TABLE]

with the definition

[TABLE]

If we assume an entrance boundary at $x_{e}$ and a reflecting boundary at $x_{c}$ for the process $x(t)\in[x_{e};x_{c}]$, then $c_{1}=0$. In this case a stationary solution of the Fokker-Planck equation (8) exists and is given by

[TABLE]

3 First passage time

We are now looking for the mean time until the process $x(t)$ reaches certain values $x_{c}$ starting at an initial value $x_{0}$. Let $\triangle\in D=(x_{l},x_{c})$ be an inner point. The first passage time of $x_{c}$ starting at $x_{0},\;\triangle<x_{0}<x_{c},$ at time $t_{0}$ is defined by $T_{x_{c}}(x_{0})=\inf\{t:x(t)=x_{c}|x(t_{0})=x_{0}\}$. The $n$-th moment $M_{n}(x_{0})$ of the first passage time $T_{x_{c}}(x_{0})$ is given by the generalized Pontryagin equation

[TABLE]

with $M_{0}(x_{0}):=1$. This equation can be solved for $n=1$ by

[TABLE]

This formula is equivalent to the formula from Karlin & Taylor (1981) on page 197.

In the following, we consider the singular case, for which the diffusion at the left boundary $x_{l}$ disappears, i.e. $\sigma^{2}(x_{l})=0$. The point $x_{l}$ is then called a singular point of the diffusion defined by the generator $\mathcal{L}$ . It is obvious that due to the singularity in equation (2) at $\sigma^{2}(x_{l})=0$, the limit behavior of equations (5) and (17) has to be determined.

For the singular case the mean first passage time $M_{1}(x_{0})$ from equation (17) can be determined. This yields our main result in the next Theorem.

Theorem 1**.**

*(Main Theorem)

Let $x_{c}$ be a regular or an exit boundary of the diffusion process $x(t)$ as defined by the Itô equation (1). Let further

(i) $x_{l}$ be an entrance boundary of the process $x(t)$

or

(ii) $S(x_{l},z]=\infty$ and $M(x_{l},z]$ $\forall x_{l}<z\leq x_{c}$.

Then the mean time until $x(t)$ reaches $x_{c}$ starting at any $x_{0}\in[x_{l},x_{c}]$ is given by*

[TABLE]

and $M_{1}(x_{0})<\infty$ $\forall\;x_{0}\in[x_{l};x_{c}]$.

Proof.

We have to show, that the limit as $\triangle\rightarrow x_{l}$ of equation (17) is given by the resulting equation (18) of Theorem 1. For the case (i) in which the singular point $x_{l}$ is an entrance boundary, we first observe that due to the definition (13) of $N_{l}(x_{l})$ we have

[TABLE]

since $N_{l}(x_{e})<\infty$. Then

[TABLE]

implies the limit

[TABLE]

Thus case (i) implies case (ii). Using conditions (ii), we can calculate the limit as $\triangle\rightarrow x_{l}$ of the first term in equation (17) by

[TABLE]

The last term of equation (17) involves the measure

[TABLE]

Therefore, that last term will not vanish in general, although $S(x_{l},x_{c}]=\infty$. The corresponding limit $Z$ is given by

[TABLE]

Because the drift $m$ and the diffusion $\sigma$ are measurable functions, we can interchange the limit with the integral and obtain

[TABLE]

Since $S(x_{e},x_{c}]=\infty$ it follows that

[TABLE]

which is the last term of equation (18).

Due to $N_{l}(x_{l})<\infty$, the process $x$ reaches the point $x_{c}$ in finite time with probability one, starting at an arbitrary $x_{0}\in[x_{l};x_{c}]$. In addition, the last term of equation (18) is also finite, because of equation (19). These facts imply $M_{1}(x_{0})<\infty$ $\forall\;x_{0}\in[x_{l};x_{c}]$. ∎

With this result, we are able to calculate the mean first passage time of nonlinear diffusion processes.

Corollary 1**.**

Let $x_{l}$ be an entrance boundary and let $x_{c}$ be a regular or an exit boundary of the diffusion process $x(t)$ as defined by the Itô equation (1). Then the mean time until $x(t)$ reaches $x_{c}$ starting at $x_{0}=x_{l}$ is finite and given by

[TABLE]

Proof.

The assertion follows from Theorem 1 because of

[TABLE]

∎

We now turn to the determination of moments of the mean first passage time $T_{x_{c}}(x_{0})$ of the process from equation (1). Starting from a point $x_{0}\in[\triangle;x_{c}]$, the $n$-th moment for reaching the regular boundary $x_{c}$ can be obtained from the generalized Pontryagin equation (16). This moment is given by

[TABLE]

where $M_{0}(x_{0}):=1$, cf. Karlin & Taylor (1981) on page 197.

For the singular case $\sigma^{2}(x_{e})=0$ at the entrance boundary $x_{e}$, we can extend Theorem 1 for the $n$-th moment of the first passage time $T_{x_{c}}(x_{0})$ as follows.

Theorem 2**.**

Let $x_{e}$ be an entrance boundary and let $x_{c}$ be a regular or an exit boundary of the diffusion process $x(t)$ as defined by the Itô equation (1). Let further $M_{n}(x_{0})$ be the $n$-th moment of the first passage time until $x(t)$ reaches $x_{c}$ starting at any $x_{0}\in[x_{e},x_{c}]$. Then

[TABLE]

where the first moment $M_{1}(x_{0})$ of the recursion is obtained from Theorem 1 and the $n$-th moment is finite, i.e. $|M_{n}(x_{0})|<\infty$ $\forall n\in\mathds{N}$.

Proof.

Let $|M_{n-1}(x_{0})|<\infty$. Then a constant $C_{n}\in\mathds{R^{+}}$ exists for any $n$, such that

[TABLE]

It follows that $|M_{n}(x_{0})|<\infty$ as well, since

[TABLE]

Because of $M_{1}(x_{0})<\infty$, we have by mathematical induction that

$M_{n}(x_{0})<\infty$ $\forall n\in\mathds{N}$. Since $M_{n}(x_{0})<\infty$, it follows that the $n$-th moment of the first passage time $T_{x_{c}}(x_{0})$ is given by equation (29). The calculations are analog to the calculations in the proof of Theorem 1. ∎

We can again consider the first passage time of $x_{c}$ starting at the entrance boundary $x_{e}$.

Corollary 2**.**

Let $x_{e}$ be an entrance boundary and let $x_{c}$ be a regular or an exit boundary of the diffusion process $x(t)$ as defined by the Itô equation (1). Then the moments of the time until $x(t)$ reaches $x_{c}$ starting at $x_{0}=x_{e}$ is finite and given by

[TABLE]

Proof.

The assertion follows from Theorem 2, since $M_{n}(x_{0})<\infty\;\forall n\in\mathds{N}$, and

[TABLE]

∎

4 Stochastic averaging

Important systems can be modeled by the perturbed Hamiltonian system

[TABLE]

with the two dimensional state space $\mathbf{Z}^{\varepsilon}:=(x,y)\in D\subset\mathds{R}^{2}$. Thereby the function $H(x,y)$ is the Hamiltonian, $d(y)$ is a damping function, $f(x,\boldsymbol{\xi}_{t})$ is a function of random excitations and $\varepsilon>0$. For the case of weakly perturbed systems of type (32) with small $\varepsilon\ll 1$, a stochastic averaging method is proposed in the following Theorem, using results by Khasminskii (1968), Borodin (1977) and Borodin & Freidlin (1995). With this method, the stochastic process of the Hamiltonian $H$ can be obtained, which is the process of total energy of the corresponding nonlinear oscillator. Such a stochastic averaging procedure was used in Dostal et al. (2012) and can be generalized as follows.

Theorem 3**.**

Let $Z^{\varepsilon}=(X,Y)\in D\subseteq\mathds{R}^{2}$ be the solution of the SDE

[TABLE]

and let the following conditions be fulfilled:

• i)

The stochastic process $\boldsymbol{\xi}_{t}:=\boldsymbol{\xi}(\omega,t)=(\xi_{1}(t),\ldots,\xi_{k}(t))^{\mathrm{T}}\in\mathds{R}^{k}$ is stationary, absolutely regular with sufficient mixing properties, and with $E\{\xi_{j}(t)\}=0$.

• i)

The functions $\mathbf{\mathrm{\mathbf{f_{0}}}}$ and $f_{1}$ satisfy certain limits in order to ensure uniqueness of the solution.

• i)

Without loss of generality, the functions $\mathbf{\mathrm{\mathbf{f_{0}}}}:D\mapsto\mathds{R}^{k}$, $f_{1}:D\mapsto\mathds{R}$ and the solution of the equation ${\dot{Y}=g(x,Y)}$ are periodic with period $T(x)$ for fixed $x$.

Let the averaging operator $\mathbb{M}$ for periodic functions $f:\mathds{R^{+}}\mapsto\mathds{R}$ with period $T$ be defined by

[TABLE]

and let $Y^{x}$ be the solution of the ordinary differential equation

[TABLE]

If the limits

[TABLE]

exist, then the process $X(\tau)$, $\tau=\varepsilon t$, converges, as $\varepsilon\rightarrow 0$, weakly on the time interval of order $\mathcal{O}(1/\varepsilon)$ to a diffusion Markov process $Z$ satisfying the Itô stochastic differential equation

[TABLE]

with the standard Wiener process $W_{\tau}$.

Proof.

The deterministic solution $Y^{x}$ of equation $\mathrm{d}Y^{x}=g(x,Y^{x})\>\mathrm{d}t$ is determined for arbitrary but fixed $x$. Then, $Y$ is replaced by $Y^{x}$ in equation (33). The assertion follows by applying the Theorem from Borodin (1977) for the resulting equation. ∎

If the functions $\mathbf{\mathrm{\mathbf{f_{0}}}}$ and $f_{1}$ are not periodic, then the procedure as described in Khasminskii (1966) has to be used. The essential result of Theorem 3 is, that the total energy $H(\textbf{Z}^{\varepsilon}(t))$ of system (32) converges in probability at a scale $\mathcal{O}(t)$, $\tau=\varepsilon t$, to the diffusion Markov process $\bar{H}(t)$ as $\varepsilon\rightarrow 0$. The resulting stochastic process is given by the Itô equation

[TABLE]

where $W_{\tau}$ is the standard Wiener process. In order to simplify the notation, we will not distinguish between the original process $H$ and the averaged process $\bar{H}$.

Remark 1**.**

In this work the function $g(X;Y)$ contained in Theorem 3 is chosen such that $g(X,Y)=\sqrt{Q}$ with $Q:=2X-2U(Y)$, whereby $U(Y):\mathds{R}\mapsto\mathds{R}$ is a continuously differentiable function. Depending on the roots of $g(X,Y)$, the drift and diffusion coefficients have to be determined piecewise for different phase space regions. More details on the dependence of the nonlinear diffusion process on $g(X,Y)$ can be found in Freidlin & Wentzell (2012).

5 Linear system with external excitation

In this section the mean and variance of first passage times is determined for amplitude crossing of a linear oscillator, which is subjected to an additive stationary wide-band random excitation. In the following sections the equations of motion of this oscillator are successively extended in order to demonstrate the developed theory. The considered linear oscillator, as also studied in Ariaratnam & Pi (1973), is given by

[TABLE]

Here, the parameter $d$ is the damping ratio, $\omega_{n}$ is the undamped natural frequency, and $\xi(t)$ is a stationary stochastic process with zero mean and spectral density

[TABLE]

where $R_{\xi\xi}(s)=E[\xi(t)\xi(t+s)]$ is the autocorrelation of the stochastic process $\xi(t)$, and $E[\cdot]$ denotes the expected value. We introduce the non-dimensional variable

[TABLE]

where $a(t)$ is the amplitude process of the oscillator (38). Then according to Ariaratnam & Pi (1973), the Itô equation for the process $r(t)$ can be written as

[TABLE]

Now the point $r=0$ is an entrance boundary, since

[TABLE]

and

[TABLE]

In equations (42) and (43), the function $Ei(\cdot)$ represents the exponential integral and $\gamma$ is Euler’s constant. We now use the Theorems 1 and 2 in order to calculate the first two moments $M_{1}(r_{0},R)$ and $M_{2}(r_{0},R)$ of the first passage time until the process $r(t)$ reaches the value $R>0$ starting at $r_{0}\in[0,R]$. From these results the mean $M(r_{0},R)=M_{1}(r_{0},R)$ and variance $V(r_{0},R)=M_{2}(r_{0},R)-(M_{1}(r_{0},R))^{2}$ of the considered first passage time are obtained. In Figure 1 the mean first passage time of the linear oscillator with additive noise Ariaratnam & Pi (1973) for reaching the boundary $R=2.2$ starting at $r_{0}$ is shown using Theorem 1 and Corollary 1 as well as the solution obtained by Ariaratnam & Pi (1973).

The corresponding variance of the first passage time obtained from Theorem 2 and from the solution provided by Ariaratnam & Pi (1973) is shown in Figure 2.

The results show, that Theorems 1 and 2 provide exact moments for the first passage time of the oscillator 38.

6 Forced and damped Mathieu oscillator

The forced and damped Mathieu oscillator is an extension of the above linear oscillator by parametric forcing. It is a basic model for many physical problems. The state space representation of the Mathieu oscillator with additive forcing $\xi_{1}(t)$ and parametric forcing $\xi_{2}(t)$ can be written as

[TABLE]

Where $t$ is time,

Rewriting equation (44) in terms of the total energy using the Hamilton function

[TABLE]

and its time derivative

[TABLE]

Combining the first equation of (44) with equation (46) and rearrange equation (45) with respect to $y$ such that

[TABLE]

we get the system

[TABLE]

An important property of equation (48) is, that the energy level $H$ changes slowly compared to the oscillations of the variable $x$. This enables the application of stochastic averaging to this system.

6.1 Mathieu oscillator excited by non-white Gaussian processes

The solution of the unforced oscillator (44) with $\varepsilon=0$, energy level $H$, and initial conditions $x(0)=0$ and $y(0)=\sqrt{2\,H}$ is given by

[TABLE]

[TABLE]

Here, $b=\sqrt{\frac{2\,H}{\alpha_{1}}}$ is the oscillation amplitude and $q=\sqrt{\alpha_{1}}$ is the natural frequency of the linear oscillator. For every energy level $H$ the oscillation period is given by

[TABLE]

Stochastic averaging of system (44) by means of theorem 3 yields the convergence of the energy $H$ for $\varepsilon\rightarrow 0$ to the solution of the one-dimensional Itô equation

[TABLE]

whereby $W_{\tau}$ is the standard Wiener process and ${\tau}=\varepsilon t$. The drift coefficient $m(H)$ in equation (52) is given by

[TABLE]

Using the trigonometric identity

[TABLE]

And the spectral density $S_{\xi_{i}\xi_{i}}(\omega)=\frac{1}{\pi}\int_{-\infty}^{0}R_{\xi_{i}\xi_{i}}(s)\cos(\omega\,s)$, we finally obtain

[TABLE]

For the diffusion coefficient $\sigma(H)$ of the equation (52) we get

[TABLE]

After evaluation of these integrals we obtain

[TABLE]

6.2 Validation of the results for the Mathieu oscillator

The results of the stochastic averaging of the Mathieu oscillator according to theorem 3 are validated using the well-known solution of the classical stochastic averaging of the stochastic linear oscillator is used. This was determined by Ariaratnam and Tam in Ariaratnam & Tam (1976). The drift and diffusion coefficients $m_{A}$ and $\sigma_{A}$ are functions of the oscillation amplitude $b$.

[TABLE]

[TABLE]

Since the drift and diffusion coefficients according to Ariaratnam & Tam (1976) are not functions of the energy $H$ but functions of the oscillation amplitude $b$, the probability densities of the oscillation amplitude $b$ are calculated in order to compare the results. The stationary solution (15) of the Fokker-Planck equation (8) using the drift and diffusion coefficients from (57) and (58) is given by

[TABLE]

where $\tilde{C}$ is a normalization constant. For the stationary probability density of the stochastic linear oscillator using the drift and diffusion coefficients (54) and (56) from the stochastic averaging of the energy by theorem 3 we obtain

[TABLE]

with a normalization constant $C$.

Using the transformation $p_{\mathrm{st}}(b)=p_{\mathrm{st}}(H)\left(\frac{\mathrm{d}}{\mathrm{d}H}b\right)^{-1}$, the probability density $p_{\mathrm{Lin}}(H)$ can be expressed a s a function of the oscillation amplitude $b$ by

[TABLE]

As can be seen from Figure 3, the two probability density functions $p_{\mathrm{Lin}}(b)$ und $p_{A}(b)$ are identical for the chosen parameters. Thus, it is obvious that the results of the two stochastic averaging methods for the stochastic linear oscillator coincide. The correspondence of these averaging methods can generally be shown by transforming the stochastic differential equation for the energy $H$ (52) into an equivalent stochastic differential equation in terms of the oscillation amplitude $b$.

7 Nonlinear system with external and parametric excitation

In this section a Duffing oscillator with softening cubic stiffness and linear, quadratic and cubic damping is investigated. Such a Duffing oscillator with the scalar state variables $x(t)$ and $y(t)$ is given by

[TABLE]

We use the Hamiltonian formalism and rewrite equation (62) in terms of the total energy defined by the Hamilton function

[TABLE]

where $\alpha_{1},\alpha_{3}>0$. Then (62) can be written as

[TABLE]

Thereby we assume $\varepsilon<<1$ to be a small parameter. With this approach we have transformed the original equation (62) to a weakly perturbed Hamiltonian system (64) excited by the stationary random processes $\xi_{1}(t)$ and $\xi_{2}(t)$. The coefficients $\nu_{1}$ and $\nu_{2}$ are additive and parametric noise intensities, respectively.

The fixed points of system (64) without dissipation and random perturbation, i.e. $\varepsilon=0$, are

[TABLE]

These saddle points $P_{1}$ and $P_{2}$ are connected by the heteroclinic orbit

[TABLE]

The considered duffing oscillator oscillates only within the region bounded by the heteroclinic orbit $\gamma$. This region is given by

[TABLE]

If the conservative system given by (62) with $\varepsilon=0$ is considered, then for every energy level $H(x,y)$ with $(x,y)\in\mathbf{D}_{\gamma}$ exactly one closed trajectory exists in the phase space of the conservative system. These trajectories correspond to the contour lines of the Hamiltonian (63) and are shown in figure 4. The time derivative of Hamiltonian for system (64) is given by

[TABLE]

Combining the first equation of (64) and equation (68), and using

[TABLE]

obtained from (63), we get the system

[TABLE]

Since as before a property of equation (70) is that the energy level $H$ changes slowly compared to the oscillations of the variable $x$, stochastic averaging can be applied.

7.1 White noise case

First we state the results for averaging system (64) subjected to white noise excitation. Therefore, let $\int_{0}^{t}\xi_{s}ds=W_{t}$, where $W_{t}$ is a standard Wiener process and $dW_{t}$ its increment. Then, using the Itô lemma, we obtain from system (64)

[TABLE]

Averaging (71) according to Khasminskii (1968) we obtain the one dimensional Itô equation

[TABLE]

where

[TABLE]

[TABLE]

Here, $T(H)$ is the period of one oscillation of the fast variable $x$ in the absence of noise and damping, i.e. $\varepsilon=0$, starting at the energy level $H$. Additionally we use

[TABLE]

For $0\leq H<\alpha_{1}^{2}/(4\alpha_{3})$ the period is given by

[TABLE]

where

[TABLE]

[TABLE]

The limits of integration $\pm b(H)$ are the points where $y=\sqrt{Q(x,H)}=0$ and the periodic orbit intersects the x-axis, i.e. $b(H)$ is the maximum value of $x$ for each energy level $H$ and is given by

[TABLE]

The function $K(k)$ is the complete elliptic integral of the first kind, cf. Byrd & Friedman (1954). The elliptic modulus $k$ is given by

[TABLE]

The integrals appearing in the equations for drift (73) and diffusion (74) exist for $0\leq H<\alpha_{1}^{2}/(4\alpha_{3})$. They can be computed in terms of complete elliptic integrals of the first and second kind, $K(k)$ and $E(k)$, respectively. Then we get

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Thus we have obtained an one-dimensional Itô equation (72) for the process of total energy $H(t)$ of system (64) subjected to white noise excitation.

7.2 Real noise case

For the non-white noise case it is necessary to determine the functions $x(t+s)$ and $y(t+s)$ in terms of $x(t)$ and $y(t)$ in order to apply stochastic averaging and obtain a closed form solution. Therefore, a solution of the differential equation (64) is needed, which can be obtained for $\varepsilon=0$ in terms of Jacobian elliptic functions. Then addition formulas for Jacobian elliptic functions can be used to eliminate the time shift and obtain the states $x$ and $y$ as functions of time $t$ only. In the case $0\leq H<\alpha_{1}^{2}/(4\alpha_{3})$ we get for $\varepsilon=0$ a solution of (64) by

[TABLE]

[TABLE]

with $b\>q=\sqrt{2H}$. The expressions for $x(\cdot)$ and $y(\cdot)$ contain the Jacobian elliptic functions $\mathrm{sn}(\cdot,k)$, $\mathrm{cn}(\cdot,k)$ and $\mathrm{dn}(\cdot,k)$, see Byrd & Friedman (1954). The elliptic modulus is the same as in the white noise case. We use the abbreviations $\mathrm{sn}:=\mathrm{sn}(qt,k),\;\mathrm{cn}:=\mathrm{cn}(qt,k),\;\mathrm{dn}:=\mathrm{dn}(qt,k),\;u:=qt.$ In addition if the subscript $s$ or $t+s$ is used, we refer to the argument $qs$ or $q(t+s)$, respectively. Applying stochastic averaging according to Theorem 3 to system (70) for $0\leq H<\frac{\alpha_{1}^{2}}{4\alpha_{3}}$, we get the one-dimensional Itô stochastic differential equation

[TABLE]

for the energy level $H$, where $\mathrm{d}W_{t}$ is a standard Wiener process. The corresponding drift $m(H)$ and diffusion $\sigma(H)$ in equation (90) are given by

[TABLE]

[TABLE]

Thereby we assume, that the involved stochastic processes $\xi_{1}(t)$ and $\xi_{2}(t)$ are stationary with mean zero and have autocorrelation function $R_{\xi_{t}\xi_{t}}(s)=E\{\xi_{t}\xi_{t+s}\}$, which approaches zero sufficiently fast as $s$ increases.

7.3 Results for first passage times

With the theory proposed in this work, the first moment of the first passage time of the Duffing oscillator with negative cubic stiffness and nonlinear damping, as given in equation (62), is determined. First, the Itô equation (90) for the averaged energy $H$ of this Duffing oscillator with equation (91) for the drift $m$ and equation (92) for the diffusion $\sigma$ is used, which is a one dimensional diffusion process. Then, Theorem 1 is applied in order to calculate the mean $M(H_{0},H_{c})=M_{1}(H_{0},H_{c})$ of the first passage time until the energy process $H(t)$ reaches the value $0<H_{c}<\frac{\alpha_{1}^{2}}{4\alpha_{3}}$ starting at $H_{0}\in(0,H_{c})$. Using the formula (79), the energy can be transformed to the oscillation amplitude $b(t)$ of the Duffing oscillator. The chosen parameter values for the softening Duffing oscillator are summarized in Table 1. These values correspond to a 200 m long ship, as described in Dostal et al. (2012), which is traveling with a velocity of 8 knots in a sea state with wave encounter angle $\chi=30^{\circ}$ a significant wave height $H_{s}=5$ m and mean wave period of $T_{1}=11.5$ s. Thereby the sea state has a JONSWAP spectral density, cf. Dostal et al. (2012). The resulting oscillation amplitude $b$ is the roll angle of the ship in degrees. The critical Energy $H_{c}=0.529$ is chosen, which corresponds to the roll oscillation amplitude $b_{c}=40$ degrees. Using the parameter values from Table 1, the functions for drift $m(b)$ and diffusion $\sigma(b)$ are shown in Figures 5 and 6, respectively. The results for the mean first passage time of the softening Duffing oscillator, for reaching the boundary $b_{c}=40$ degrees starting at $b_{0}$ are shown in Figure 7 using Theorem 1 and Corollary 1. The computation time of the corresponding formula (18) from Theorem 1 using standard numerical quadrature is about 15 seconds on an ordinary desktop computer. The mean first passage time obtained from Corollary 1 is only correct, if the starting oscillation amplitude $b_{0}$ is very close to zero. For higher values of the starting oscillation amplitude $b_{0}$, the mean first passage times have to be determined using Theorem 1.

8 Conclusions

The moments of the first passage time are obtained for a one-dimensional nonlinear diffusion processes with an entrance boundary by examining the boundary behavior. Thereby previous theorems were extended, such that they are valid for arbitrary initial points of the considered diffusion. The mean and variance of the first passage time to reach the boundary of a domain are validated with known analytical formulas, which perfectly match. Results of the first passage times for a softening Duffing oscillator are obtained as well, which are important for the determination of dangerous ship roll dynamics in ocean waves. This shows, that the proposed theory is applicable for important problems. The necessary computation time is very low, since the proposed analytical expressions for the moments of the first passage times only involve integrals, which can be evaluated using standard quadrature formulas.