A Feynman-kac Formula Approach for Computing Expectations and Threshold Crossing Probabilities of Non-smooth Stochastic Dynamical Systems

TL;DR

This paper introduces a Feynman-Kac formula-based numerical method for efficiently computing expectations and crossing probabilities in non-smooth stochastic systems like elasto-plastic and obstacle problems, avoiding extensive simulations.

Contribution

It develops a novel PDE approach for non-smooth stochastic systems with complex boundary conditions, providing a computational alternative to probabilistic simulations.

Findings

Effective numerical solutions for variational inequality PDEs

Accurate computation of threshold crossing probabilities

Application to engineering mechanics problems

Abstract

We present a computational alternative to probabilistic simulations for non-smooth stochastic dynamical systems that are prevalent in engineering mechanics. As examples, we target (1) stochastic elasto-plastic problems, which involve transitions between elastic and plastic states, and (2) obstacle problems with noise, which involve discrete impulses due to collisions with an obstacle. We formally introduce a class of partial differential equations related to the Feynman-Kac formula, where the underlying stochastic processes satisfy variational inequalities modelling elasto-plastic and obstacle oscillators. We then focus on solving them numerically. The main challenge in solving these equations is the non-standard boundary conditions which describe the behavior of the underlying process on the boundary. We illustrate how to use our approach to compute expectations and other statistical quantities, such as the asymptotic growth rate of variance in asymptotic formulae for threshold crossing probabilities.

Tables

Table 1. Table 1 . Results using the MC approach with different T 𝑇 T and δ ​ t 𝛿 𝑡 \delta t for the approximation of the probability ( E 1 subscript 𝐸 1 E_{1} ) (left) and the growth rate of the plastic deformation ( E 3 subscript 𝐸 3 E_{3} ) (right). Here P Y = 0.25 subscript 𝑃 𝑌 0.25 P_{Y}=0.25 .

$T$	$\delta t={10}^{-2}$	$\delta t={10}^{-3}$	$\delta t={10}^{-4}$
2	0.518	0.487	0.474
4	0.519	0.492	0.478
8	0.519	0.489	0.477
16	0.520	0.487	0.481
32	0.518	0.491	0.478
64	0.519	0.490	0.476

Table 2. Table 2 . Results using PDE and MC method approaches for the elasto-plastic problem ( 3 ) for different plastic yields P Y subscript 𝑃 𝑌 P_{Y} : Probability ( E 1 subscript 𝐸 1 E_{1} ) of plastic state for T → ∞ , ( T ≥ 5 ) → 𝑇 𝑇 5 T\to\infty,(T\geq 5) . Asymptotic growth rate ( E 3 subscript 𝐸 3 E_{3} ) of the plastic/total deformation for T → ∞ , ( T ≥ 150 ) → 𝑇 𝑇 150 T\to\infty,(T\geq 150) .

	$P_Y$	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1.0
($E_1$)	PDE	0.638	0.521	0.430	0.354	0.289	0.234	0.187	0.148	0.115	0.088
$T\to \infty$	MC	0.639	0.521	0.429	0.352	0.286	0.230	0.184	0.144	0.112	0.085
($E_3$)	PDE	0.776	0.589	0.440	0.325	0.238	0.173	0.125	0.089	0.064	0.045
$T\to \infty$	MC	0.799	0.607	0.457	0.337	0.248	0.182	0.132	0.096	0.070	0.050

Table 3. Table 3 . Results using PDE and MC method approaches for the obstacle problem ( 5 ) for different position of the obstacle P O subscript 𝑃 𝑂 P_{O} : Mean kinetic energy ( E 2 ′ superscript subscript 𝐸 2 ′ E_{2}^{\prime} ) for T → ∞ , ( T ≥ 5 ) → 𝑇 𝑇 5 T\to\infty,(T\geq 5) . Asymptotic growth rate ( E 3 subscript 𝐸 3 E_{3} ) of the integral of the deformation for T → ∞ , ( T ≥ 150 ) → 𝑇 𝑇 150 T\to\infty,(T\geq 150) .

	$P_Y$	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1.0
($E_2^{\prime }$)	PDE	0.179	0.250	0.298	0.347	0.364	0.388	0.409	0.426	0.441	0.453
$T\to \infty$	MC	0.179	0.250	0.297	0.334	0.364	0.389	0.408	0.423	0.437	0.450
($E_3^{\prime }$)	PDE	0.00148	0.0097	0.028	0.062	0.103	0.158	0.223	0.294	0.370	0.447
$T\to \infty$	MC	0.00149	0.0097	0.028	0.060	0.103	0.158	0.222	0.295	0.370	0.447

Equations

$$\dot{X}_{t}=F(X_{t},Y_{t}),\quad\dot{Y}_{t}=G(X_{t},Y_{t})+{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\mathbf{H}_{t}}}+\dot{W}_{t},\quad\forall t>0$$

$$A\triangleq \mathbb{E}\left(f(X_T,Y_T)+{\int }_0^{T}g(X_{\tau },Y_{\tau })d\tau ,\right),B\triangleq \mathbb{E}\left({\int }_0^{\infty }{e}^{-\lambda \tau }g(X_{\tau },Y_{\tau })d\tau \right),$$

$$C\triangleq\lim_{T\to\infty}\mathbb{E}\left(f(X_{T},Y_{T})\right),\quad{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}D\triangleq\mathbb{P}\left(\max_{0\leq t\leq T}\left|\int_{0}^{t}f(X_{\tau},Y_{\tau})\hbox{d}\tau\right|\geq b\right),}}$$

$$A^{\prime }$$

$$B^{\prime }$$

$$C^{\prime }$$

$$F(x,y)\triangleq \frac{\partial \mathcal{H}}{\partial x}(x,y),G(x,y)\triangleq \frac{\partial \mathcal{H}}{\partial y}(x,y)+\eta (x,y)\frac{\partial \mathcal{H}}{\partial x}(x,y),(x,y)\in {\mathbb{R}}^2,$$

$${\dot{Z}}_t$$

$${\dot{\Delta }}_t$$

$${{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\mathbf{H}_{t}}}=kZ_{t},\quad k>0.$$

$${\Delta }_t={\int }_0^t{1}_{\{|Z_s|=P_Y\}}{Y}_{s}ds.$$

$$\forall t,\forall |\varphi |\le P_Y,d{Y}_t=-(c_0{Y}_t+k{Z}_t)dt+d{W}_t,(d{Z}_t-Y_{t}dt)(\varphi -Z_t)\ge 0$$

$$E_1$$

$$E_2$$

$$E_3$$

$$\mathbb{E}f(Z_T,Y_T)\phi (Z_{T+h},Y_{T+h}),$$

$$E_4$$

$$\underset{T\to \infty }{\lim }\mathbb{P}\left(|Z_{T+h}|=P_Y,|Z_T|=P_Y\right),\underset{T\to \infty }{\lim }\mathbb{E}{Y}_T^2{Y}_{T+h}^2,$$

$$E_5$$

$${{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\mathbf{H}_{t}=\dot{I}_{t}.}}$$

$$I_t\triangleq \sum _{0\le s\le t}\left({\dot{X}}_{s+}-{\dot{X}}_{s-}\right)1_{\{|X_s|=P_O\}}.$$

$$X_t=x+{\int }_0^t{Y}_{s}ds,Y_t=y-{\int }_0^t(c_0{Y}_s+k{X}_s)ds+W_t-(1+e)\sum _{0\le s\le t}{\dot{Y}}_{s-}{1}_{\{|X_s|=P_O\}}$$

$$d{X}^0=Y^{0}dt,d{Y}^0=-(c_0{Y}^0+k{X}^0)dt+dW$$

$$E_1^{\prime }$$

$$E_2^{\prime }$$

$$E_3^{\prime }$$

$$E_4^{\prime }$$

$$\underset{T\to \infty }{\lim }\mathbb{P}\left(f(X_{T+h},Y_{T+h})\le ϵ,f(X_T,Y_T)\le ϵ\right),\underset{T\to \infty }{\lim }\mathbb{E}{Y}_T^2{Y}_{T+h}^2,$$

$$E_5^{\prime }$$

$$\{\begin{array}{ll}& \text{d}{X}_t^{ϵ}={ϵ}^{2}f(X_t^{ϵ},Y_t^{ϵ})\text{d}t+ϵg(X_t^{ϵ},Y_t^{ϵ})\text{d}{W}_t,X_0^{ϵ}=x,\\ & \text{d}{Y}_t^{ϵ}=h(X_t^{ϵ},Y_t^{ϵ})\text{d}t+\phi (X_t^{ϵ},Y_t^{ϵ})\text{d}{W}_t,Y_0^{ϵ}=y\end{array}$$

$$\text{d}{\bar{X}}_t={ϵ}^2\bar{f}({\bar{X}}_t)\text{d}t+ϵ\bar{g}({\bar{X}}_t)\text{d}{W}_t,{\bar{X}}_0^{ϵ}=x,$$

$$dY(t)=h(x,Y_t)\text{d}t+\phi (x,Y_t)\text{d}{W}_t,Y_0^{ϵ}=y.$$

Full text

Laurent Mertz***ECNU-NYU Institute of Mathematical Sciences, NYU Shanghai, Shanghai, China, Georg Stadler†††Courant Institute of Mathematical Sciences, New York University, New York, USA, Jonathan Wylie‡‡‡Department of Mathematics, City University of Hong Kong, Hong Kong, China

1. Motivations and Goal

In this paper, we propose a computational alternative to probabilistic simulations for a certain type of non-smooth stochastic dynamical systems, namely

[TABLE]

with initial condition $(X_{0},Y_{0})=(x,y)\in\mathbb{R}^{2}$. Here, $(X_{t},Y_{t})\in\mathbb{R}^{2}$ is the state variable at time $t$, dots denote derivatives with respect to time, $F(x,y)$ and $G(x,y)$ are deterministic functions, $\dot{W}$ represents a white noise random forcing in the sense that $W$ is a Wiener process, and $\mathbf{H}$ is a functional depending on interactions of the state variable at boundaries and interfaces. We are interested in statistical quantities that characterize and predict the behavior of $\{(X_{t},Y_{t}),t\geq 0\}$ such as means, moments, correlations and probabilities of crossing thresholds. These are of the form

[TABLE]

[TABLE]

and

[TABLE]

Here and in the remainder of the paper, $\mu,\lambda$ are positive numbers, $b$ is a given threshold, $T>0$ is a given time, $h\geq 0$ and $f,g,\varphi,\psi$ are continuous functions. For cases in which $\mathbf{H}\equiv 0$ and $(F,G)$ satisfy appropriate smoothness conditions, a natural setting for characterizing such quantities with partial differential equations (PDEs) is to use the Feynman-Kac formula (FKf) [45]. Furthermore, if in addition $F$ and $G$ can be written in terms of a Hamiltonian structure $\mathcal{H}(x,y)$, in the sense

[TABLE]

where $\eta$ is a well behaved function, then the process $(X_{t},Y_{t})$ belongs to the class of Stochastic Hamiltonian Systems (SHS) [48] for which quantities of type $C$ and $C^{\prime}$ are well defined. Applications that make use of this framework abound in many areas of science and engineering, e.g., finance, chemistry, biology, neuroscience, economy and mechanics.

There are many important applications involving non-smooth dynamics in which quantities of the form A,B,C,{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}D}} and $A^{\prime},B^{\prime},C^{\prime}$ are of interest. We have in mind problems involving interactions with boundaries, constraints, phase transitions or hysteresis. An important class of examples includes elasto-plasticity (EP) problems with random forcing [26, 20] where the dynamics takes into account phase transitions between elastic and plastic states. A second class of examples includes stochastically driven obstacle problems [1, 18] where the dynamics has to take into account instantaneous collisions with an obstacle. The constitutive models are shown in Figure 1. A number of authors have applied innovative techniques to determine statistics of non-smooth oscillators [44, 19, 50, 49, 51, 46, 36, 34, 29, 31]. They used methods such as stochastic averaging, statistical linearization and numerical path integration. However, some of these methods are based on approximations of the underlying mechanical system, whereas the methodology that we present in this paper does not make any such approximations and therefore provides a rigorous basis for calculating the above-mentioned quantities. A concise review on the stochastic engineering dynamics literature is provided in section 2.

From a mathematical viewpoint, these problems can be reformulated in terms of (degenerate) stochastic variational inequalities (SVIs) [10, 3]. With this formulation, the evolution of the state variables is Markovian and Kolmogorov Equations (KEs) can be derived.

Except for very few non-smooth dynamical systems, where analytical expressions are available, in general one has to resort to computational techniques to determine A,B,C,{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}D}} and $A^{\prime},B^{\prime},C^{\prime}$. The most straightforward computational methods are direct probabilistic simulations. Such techniques are simple to implement and widely used. For low-dimensional state variables, they are less efficient than techniques based on partial differential equations (PDEs), provided the latter are available. For the two non-smooth stochastic processes targeted in this paper, the presence of constraints or phase transitions (governed by a variational inequality structure) leads to non-standard boundary conditions in the PDEs for FKf. These boundary conditions characterize the behavior of the underlying process on the boundary. This paper is concerned with the numerical treatment of such non-standard PDEs and their application to compute the quantities of A,B,C{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0},D}} and $A^{\prime},B^{\prime},C^{\prime}$.

We point out that the idea of solving non-standard PDEs for FKf in the context of a SVI has been used before [5]. However, the technique in [5] is specific to the elasto-perfectly-plastic oscillator excited by white noise, for which the non-standard condition has to be satisfied in a finite (namely two) number of points. It has no natural extension to problems of obstacle type because for these problems the non-standard condition must be satisfied on a continuous set of points.

In the remainder of this section, we give the definitions of the elasto-plastic and obstacle models in presence of noise, and we discuss the quantities we are interested in.

1.1. Elasto-plastic problem with noise

A basic prototype for modeling mechanical structures that admit permanent deformation under vibrations is the elastic-plastic oscillator [26]. The dynamics focuses on two quantities: the total deformation, $X_{t}$, supported by the structure when subjected to vibrations, and its velocity, $\dot{X}_{t}$. For general elasto-plastic problems, the nonlinear functional $\mathbf{H}$ in (1) is a restoring force arising from the structure. The exact form of $\mathbf{H}$ depends on the particular structure in question. The nonlinearity in such models comes from the switching of regimes from an elastic phase to a plastic one, or vice versa. For the elasto-perfectly-plastic oscillator (EPPO) model, the irreversible (plastic) deformation $\Delta$ and the reversible (elastic) deformation $Z$ at time $t$ satisfy

[TABLE]

where $X_{t}=Z_{t}+\Delta_{t}$. Typically, $|Z_{t}|$ is bounded by a given threshold $P_{Y}$ at all times $t$, the system is in the plastic phase when $|Z_{t}|=P_{Y}$ and the elastic phase when $|Z_{t}|<P_{Y}$. Here, $P_{Y}$ is an elasto-plastic bound, known as the “Plastic Yield” in the engineering literature. We assume the force $\mathbf{H}$ is a linear function of $Z$ (while $\Delta$ remains constant) in the elastic phase and a constant (while $Z$ remains constant at $\pm P_{Y}$) in the plastic phase as follows:

[TABLE]

The permanent (plastic) deformation at time $t$ can then be written as

[TABLE]

We consider the case of linearly damped spring for which $F(x,y)=y$ and $G(x,y)=-c_{0}y$.

SVI framework

It has been shown that the dynamics of such a nonlinear oscillator can be described mathematically by means of SVIs [10]. The dynamics is then described by the pair $(Y_{t},Z_{t})$ that satisfies

[TABLE]

and appropriate initial conditions for $Y_{0}$ and $Z_{0}$ must be prescribed. Here, whereas $Y_{t}=\dot{X}_{t}$, $X_{t}$ is not involved in the dynamics.

Characterizing the statistics of the state variable using FKf is interesting for engineering purposes. To illustrate the robustness and efficiency of our approach we will compute five important quantities in this paper:

• •

the probability of finding the system in the plastic state

[TABLE]

• •

the mean kinetic energy

[TABLE]

• •

the variance of both the plastic deformation and the total deformation

[TABLE]

• •

any correlation structure between $(Y_{T},Z_{T})$ and $(Y_{T+h},Z_{T+h})$ of the form

[TABLE]

where here we will target

[TABLE]

and

[TABLE]

• •

probabilities of threshold crossings for the total deformation

[TABLE]

The probability of finding the system in the plastic state and the mean kinetic energy, at a given time $T$, are quantities of type $A$ whereas, at large time, they are of type $C$. Also, their Laplace transforms belong to type $B$. The variance of the plastic deformation and of the total deformation, at a given time $T$, are quantities of type $A^{\prime}$ and their asymptotic growth belongs to the type $C^{\prime}$. Any correlation structure of $(Y,Z)$ between two different instants $T$ and $T+h$ is of type $A^{\prime}$ and its corresponding double Laplace transform (with parameters $\lambda$ and $\mu$) is of type $B^{\prime}$. The probabilities of threshold crossings for the plastic deformation are of type $D$. These quantities are challenging to compute because analytic formulas are not available and probabilistic simulations are computationally expensive.

1.2. Obstacle problem with noise

One of the simplest models exhibiting a vibro-impact motion can be expressed as a single degree of freedom oscillator constrained by obstacles [1] located at $|X|=P_{O}$ (position of the obstacle). It is common in the engineering literature to formulate the dynamics of a stochastic obstacle oscillator in terms of a stochastic process $X_{t}$, the oscillator displacement. For general obstacle problems, the nonlinearity in such models comes from the collisions. If at a time $t$, the state hits the obstacle with incoming velocity $\dot{X}_{t}$, it immediately bounces back with velocity $-e\dot{X}_{t}$, that is, $\dot{X}_{t+}=-e\dot{X}_{t-}$, $e\in[0,1]$. The nonlinear functional $\mathbf{H}$ in (1) is a force keeping track of the past discrete impulses due to collisions with an obstacle, here

[TABLE]

The net impulse process, which keeps track of the sum of all past impulses, at time $t$ can be written as

[TABLE]

For simplicity, we consider the case of a linearly damped spring: $F(x,y)=y$ and $G(x,y)=-c_{0}y-kx$.

Reflected Langevin process framework

From a mathematical viewpoint (using a stochastic differential framework), the dynamics of such a nonlinear oscillator can be described in the framework of a Reflected Langevin process [23, 24] of the form

[TABLE]

where $Y_{t}\triangleq\dot{X}_{t}$.

Remark 1**.**

*When $e=1$, Equation (5) can be formulated in the framework of SVIs [3] as follows $\forall t,\>|X_{t}|\leq P_{O},\>\forall|\varphi|\leq P_{O},\>(\hbox{d}Y_{t}+(c_{0}Y_{t}+kX_{t})\hbox{d}t-\hbox{d}W_{t})(\varphi-X_{t})\geq 0.$ *

Comments

The implementation of the impact rule can be explained using stopping times. We start by defining $\tau_{0}\triangleq 0$ and $\{(X_{t}^{0},Y_{t}^{0}),t\geq 0\}$ to be the solution of the unconstrained problem

[TABLE]

and define $\tau_{1}\triangleq\inf\{t>\tau_{0},\>|X_{t}^{0}|=P_{O}\}$. For $t\geq\tau_{1}$, we define $\{(X_{t}^{1},Y_{t}^{1}),t\geq\tau_{1}\}$ to be the solution of (6) with the initial condition $X_{\tau_{1}}^{1}\triangleq X_{\tau_{1}}^{0}$ and $Y_{\tau_{1}}^{1}\triangleq-eY_{\tau_{1}}^{0}$ and then we define $\tau_{2}\triangleq\inf\{t>\tau_{1},\>|X_{t}^{1}|=P_{O}\}$. Knowing $\tau_{n},X^{n},Y^{n}$, we can recursively define $\tau_{n+1},X^{n+1},Y^{n+1}$. Thus, the whole process can be defined on each interval $[\tau_{n},\tau_{n+1})$ by setting $(X,Y)\triangleq(X^{n},Y^{n})$. It should be pointed out that the presentation remains formal as there are mathematical subtleties behind the definition of these stopping times. An existence and uniqueness result in a weak sense (see for instance, definition 3.1 page 300 in [25]) for an equation similar to (5) can be found in the works of Jacob [23, 24]. To be more precise, it concerns an integrated Wiener process constrained to stay in $[0,\infty)$ by a partially elastic boundary at [math]. As explained in Jacob’s works, even if the aforementioned process enjoys the property of local pathwise existence and uniqueness away from the obstacle, global pathwise uniqueness results do not hold. Indeed, even for the deterministic problem it has been shown in the work of Ballard [2] that the only case in which there exists a unique solution is when the forcing is an analytic function (locally given by a convergent power series), and this condition cannot be relaxed.

For the obstacle problem, we consider the following five quantities:

• •

the probability of finding the system in the neighborhood of the obstacle with a low velocity,

[TABLE]

where $f(x,y)\triangleq\sqrt{(|x|-P_{O})^{2}+y^{2}}$.

• •

the mean kinetic energy

[TABLE]

• •

the variance of the integral of the displacement

[TABLE]

Up to a multiplicative constant, this is equivalent to computing the variance of the change in momentum due to the restoring force $kX$.

• •

any correlation structure between $(X_{T},Y_{T})$ and $(X_{T+h},Y_{T+h})$ of the form

[TABLE]

and

[TABLE]

• •

probabilities of threshold crossings for the integral of the displacement

[TABLE]

An explanation similar to what was provided above for the elasto-plastic problem (in terms of plastic state, variance of the plastic deformation and correlations) applies to the obstacle problem (in terms of mean kinetic energy, variance of the integral of the displacement, correlations and probabilities of threshold crossings for the integral of the displacement).

1.3. Approach and overview

The PDEs related to $A,B,C$ and $A^{\prime},B^{\prime},C^{\prime}$ including quantities $(E_{i})$, $(E_{i}^{\prime})$, $1\leq i\leq 4$ for the elasto-plastic and obstacle problems are non-standard boundary value problems. For this type of non-standard PDEs, (for instance see Section 3.2), only partial existence and uniqueness results are available, mainly for the case of the EPPO problem with noise [10, 11, 9, 6]. This is because standard PDE theory techniques do not apply due to the non-standard boundary conditions and the degeneracy of these problems. Therefore, in this essay, we mostly study the behavior of these PDEs numerically in order to gather insight in the solution behavior and we compare with probabilistic approximations in order to conjecture whether solutions exist or not. For this purpose, we solve the PDE problems on sequences of grids with increasing resolution and monitor the behavior of the numerical solutions. If the numerical solutions computed from differently accurate discretization yield a converging behavior as the mesh is refined, we can conjecture that the continuous problem has a solution.

Monte Carlo methods are used to compare the PDE results with probabilistic results. These comparisons increase our confidence in the solution of the PDE problem, allow us to study approximation errors in both schemes, and establish a direct connection between the SVIs and the PDE problems.

1.4. Organization of the paper

A review of established stochastic engineering methodologies is given in section 2. In section 3, we show the connection between the quantities $A$ and $A^{\prime}$ to the solution of non-standard parabolic problems, the functions $B$ and $B^{\prime}$ to the solution of non-standard elliptic problems and the functions $C$ and $C^{\prime}$ to the solution of non-standard Poisson problems. In Section 4, an asymptotic formula for D is presented when $T$ and $b$ are large enough. This formula relies on quantities of type $C^{\prime}$. In Section 5, we present a numerical approach for solving non-standard PDEs. The method is first presented and applied to the elasto-plastic problem. Then, it is applied to the obstacle problem. Numerical results are presented. In Section 6, we compare the approach proposed in this paper with previous techniques employed for a white noise EPPO. Finally, in Section 7, broader impacts of the present method are discussed for promising engineering applications.

2. Review of stochastic engineering methodologies

We briefly review three main established techniques, namely i) stochastic averaging, ii) statistical linearization and iii) numerical path integration. The notation used in this section only applies to this section and not beyond.

2.1. Statistical linearization

The statistical linearization (SL) method, proposed by [15], has proved to be a very useful approximation technique over the years. See [12] for the mathematical validity of the SL method in some cases. In the context of a basic illustrative example, it consists in replacing an equation of the form $g(Y_{t},\dot{Y}_{t},\ddot{Y}_{t})=X_{t}$ with an equivalent form $m^{\star}\ddot{Y}_{t}+c^{\star}\dot{Y}_{t}+k^{\star}Y_{t}=X_{t}$ where $Y,\dot{Y},\ddot{Y}$ are, respectively, the displacement response, the velocity and the acceleration, $g$ is a (possibly non-linear) function and $X$ is an input forcing. The parameters (that can depend on the time) $(m^{\star},c^{\star},k^{\star})$ are determined by minimizing the error $\epsilon(c,m,k)\triangleq g(Y_{t},\dot{Y}_{t},\ddot{Y}_{t})-m\ddot{Y}_{t}-c\dot{Y}_{t}-kY_{t}$ as follows: $\|\epsilon(c^{\star},m^{\star},k^{\star})\|=\min_{c,m,k\in\mathbb{R}}\|\epsilon(c,m,k)\|$ where $\|.\|$ is a convenient norm. Further details can be found in [36], with an introduction of SL for simple systems in Chapter 5 and analysis of systems with multiple degrees of freedom in Chapter 6. The method can be employed to deal with complex systems (including hysteretic elements) having many degrees of freedom and a broad class of excitation (even non-stationary), SL usually gives reasonably good results when non-linear effects are present.

2.2. Stochastic averaging

The method of stochastic averaging, originally introduced by [47], provides approximate solutions to problems involving the vibration response of lightly damped systems to broad-band random excitation. See [27, 28, 39, 40] for mathematical rigorous foundations. For illustration, we consider one-dimensional variables and the presentation we use is inspired by [16]. Consider a system of the form

[TABLE]

for $0<\epsilon\ll 1$ and $f,g,h,\varphi$ functions with appropriate conditions. Here, $W$ is a Wiener process. The variables $X^{\epsilon}_{t}$ and $Y^{\epsilon}_{t}$ are, respectively, called slow and fast components. Engineers are interested in the behavior of $X^{\epsilon}_{t}$ on intervals of order $\epsilon^{-1}$ since on those time scales the most significant changes occur. The averaging principle states that a good approximation of the slow component can be obtained by averaging over the fast components in the sense that the trajectory $X^{\epsilon}_{t}$ is replaced by the solution $\bar{X}$ of

[TABLE]

where $\bar{f}(x)=\lim\limits_{T\to\infty}\frac{1}{T}\mathbb{E}\int_{0}^{T}f(x,Y_{t})\textup{d}t$ and $\bar{g}^{2}(x)=\lim\limits_{T\to\infty}\frac{1}{T}\mathbb{E}\int_{0}^{T}g^{2}(x,Y_{t})\textup{d}t$ and

[TABLE]

In the dynamics of $Y_{t}$, $x$ is a frozen parameter. The stochastic averaging method is particularly useful in dealing with situations with fairly light non-linear damping. Details on the application of this method to systems of the form $\ddot{X}+\epsilon^{2}\psi(X,\dot{X})+kX=\epsilon\dot{W}$ are presented in [37]. To apply the averaging principle the joint response $(X,\dot{X})$ needs to be transformed in $X_{t}=a_{t}\cos(\omega_{0}t+\varphi_{t})$ and $\dot{X}_{t}=-a_{t}\omega\cos(\omega_{0}t+\varphi_{t})$. This leads to a stochastic differential equation for $(a_{t},\varphi_{t})$ which can be seen as the slow component in the equation above. Details can be found in Section 3 of [37]. The stochastic averaging method can cope with systems having hysteretic features such as the elasto-plastic and impact problems.

2.3. Numerical path integration

The path integration (PI) method, introduced in [34], is a step-by-step calculation of the joint probability density function (PDF) of a set of state space variables describing a white noise excited nonlinear dynamical system. Precisely, the PDF is computed at a given time by applying the Chapman-Kolmogorov equation when the response PDF and the transition PDF are known at an earlier time. The PI method is an efficient approximation for solving the Fokker-Planck equation and for providing the stationary response of the underlying dynamical system. The method exploits short time transition probability density functions which can be sharply approximated by explicit formulae (based on a local Gaussian behavior) even for hysteretic systems such as elasto-plastic or impact problems.

2.4. Comments

The statistical linearization and stochastic averaging methods rely on an approximation of the underlying stochastic process, where in our method no such approximation is needed. However numerical path integration methods deal with the original process. Our approach relies heavily on the PDEs related to the Feynman-Kac (FK) formula and the backward Kolmogorov equations whereas the PI method uses the Chapman-Kolmogorov equation which ultimately lead to the Forward Kolmogorov equation. Hence PI and FK methods can be seen as being the dual of each other. In contrast with the PI method, our method extends naturally (by adapting the infinitesimal generator) to solving a broad range of problems. These include optimal stopping and stochastic optimal control problems via free boundary value problem and HJB equations with non-standard boundary condition, see [32]. Our approach can be combined with existing methodologies to investigate a number of practically important quantities. For example, following the steps of [43] and [14] it can be used to study the power spectrum density for stationary processes of the form

[TABLE]

where $X,Y$ is the response variable of one of the systems studied in this paper (elasto-plastic or impact system) and the notation $|.|$ stands for the modulus of a complex number. Some details are given in the conclusions and perspectives section.

3. The Partial Differential Equations for $A,B,C$ and $A^{\prime},B^{\prime},C^{\prime}$

It is possible to relate the functions $A$ and $A^{\prime}$ to the solution of parabolic problems, the functions $B$ and $B^{\prime}$ to the solution of elliptic problems and the functions $C$ and $C^{\prime}$ to the solution of Poisson problems. In the first part of this section, we present these problems in the case of $\mathbf{H}=0$. Then, in the second part, we provide a formal presentation of the corresponding problems for the two non-smooth problems targeted in this paper.

Notation

For $T>0$ and a domain $\Omega$ of $\mathbb{R}^{2}$, we use the notation $C^{\star}(\Omega\times[0,T])$ for the set of continuous functions on $\Omega\times[0,T]$ that are $C^{1}$-regular with respect to $x$, $C^{2}$-regular with respect to $y$ and $C^{1}$-regular with respect to $t$. We use the notation $C^{\star}(\Omega)$ for the set of continuous functions on $\Omega$ that are $C^{1}$-regular with respect to $x$ and $C^{2}$-regular with respect to $y$. We use the generic notation $\Gamma$ for

[TABLE]

because it helps to write the quantities $A,B,C$ and $A^{\prime},B^{\prime},C^{\prime}$ in a compact form.

We use the following notations below

[TABLE]

[TABLE]

[TABLE]

[TABLE]

3.1. The case $\mathbf{H}=0$

The proofs of connecting the PDEs below and stochastic processes can be found in the Appendix for the reader’s convenience. We use the notation

[TABLE]

Backward-in-time parabolic problems

Consider functions $u\in C^{\star}(\mathbb{R}^{2}\times[0,T]),v\in C^{\star}(\mathbb{R}^{2}\times[0,T+h]),w\in C^{\star}(\mathbb{R}^{2}\times[0,T])$. Let $\lambda,\mu\geq 0.$ Assume that $u,v,w$ satisfy $\mathcal{A}_{0}(\lambda,u)$, $\mathcal{A}_{0}(\mu,v)$, $\mathcal{A}_{0}(\lambda+\mu,w)$, respectively, and that

[TABLE]

and

[TABLE]

Then we have

[TABLE]

As a corollary,

[TABLE]

where $w$ solves the problem (8) with $v=u$ and $h=0$.

Degenerate elliptic problems

Let $\lambda,\mu\geq 0$ be two numbers and $g,\psi$ be two functions such that

[TABLE]

Consider functions $u_{\lambda},v_{\mu},w_{\lambda+\mu}\in C^{\star}(\mathbb{R}^{2})$. Assume that $u_{\lambda},v_{\mu},w_{\lambda+\mu}$ satisfy technical conditions $\mathcal{A}_{1}(\lambda,u_{\lambda}),\mathcal{A}_{1}(\mu,v_{\mu})$, $\mathcal{A}_{1}(\lambda+\mu,w_{\lambda+\mu})$, $\mathcal{A}_{1}(\lambda,\mu,\frac{\partial u_{\lambda}}{\partial y}\frac{\partial v_{\mu}}{\partial y})$, $\mathcal{A}_{1}(\lambda+\mu,u_{\lambda}v_{\mu})$ $\mathcal{A}_{3}(\lambda,\mu,u_{\lambda},\psi)$, $\mathcal{A}_{3}(\mu,\lambda,v_{\mu},g)$ and

[TABLE]

Then we have

[TABLE]

Long time behavior and degenerate Poisson problems

Assume that $(X_{t},Y_{t})$ has a unique invariant probability measure $\nu$ on $\mathbb{R}^{2}$ in the sense that for any continuous function $f$ satisfying $\nu(|f|)<\infty$, we have $\mathbb{E}f(X_{t},Y_{t})=\nu(f)$ provided that $(X_{0},Y_{0})$ is distributed according to $\nu$.

Consider functions $U,V\in C^{\star}(\mathbb{R}^{2})$. Assume that $g$ and $\psi$ satisfy $\nu(|g|)<\infty,\nu(|\psi|)<\infty$ and $U,V$ satisfy

[TABLE]

together with

[TABLE]

Then we have

[TABLE]

3.2. Non-standard PDEs of SVIs: elasto-plastic and obstacle problems

In a similar manner to the results above for the case $\mathbf{H}=0$, here we present non-standard PDEs related to SVIs modeling the elasto-plastic (2) and obstacle (4) problems. First, for each problem, we give a description of the infinitesimal generator of the corresponding process.Then, we only present the non-standard backward-in-time parabolic problems in analogy to Section 3.1. For the non-standard elliptic degenerate and Poisson problems, the idea remains basically the same as what is presented above.

Without loss of generality, we can assume here that $P_{Y}=P_{O}=1$, define

[TABLE]

and

[TABLE]

In the elasto-plastic problem, using Ito’s lemma, the generator of $(Z_{t},Y_{t})$ is defined on any function $\phi\in C^{\star}(\bar{D})$ and satisfies

[TABLE]

In the obstacle problem, (this is formal) the generator of $(X_{t},Y_{t})$ is defined on any function $\phi\in C^{\star}(\bar{D})$ such that

[TABLE]

satisfies

[TABLE]

Non-standard problems for the elasto-plastic problem

Consider functions $u\in C^{\star}(D\times[0,T]),v\in C^{\star}(D\times[0,T+h]),w\in C^{\star}(D\times[0,T])$. Assume that $u,v,w$ satisfy $\mathcal{A}_{0}(\lambda,u)$, $\mathcal{A}_{0}(\mu,v)$, $\mathcal{A}_{0}(\lambda+\mu,w)$, respectively, (w.r.t to the solution of (3))

[TABLE]

[TABLE]

[TABLE]

then the corresponding $A$ and $A^{\prime}$ satisfy

[TABLE]

Non-standard problems for the obstacle problem

Consider functions $u\in C^{\star}(D\times[0,T]),v\in C^{\star}(D\times[0,T+h]),w\in C^{\star}(D\times[0,T])$. Assume that $u,v,w$ satisfy $\mathcal{A}_{0}(\lambda,u)$, $\mathcal{A}_{0}(\mu,v)$, $\mathcal{A}_{0}(\lambda+\mu,w)$, respectively, (w.r.t to the solution of (5))

[TABLE]

[TABLE]

[TABLE]

then the corresponding $A$ and $A^{\prime}$ satisfy

[TABLE]

4. Asymptotic formula for D

For every function $f$ on $\mathbb{R}^{2}$, define $\Delta_{t}^{EP}(f)\triangleq\int_{0}^{t}f(Z_{s},Y_{s})\hbox{d}s\>\mbox{ and }\>\Delta_{t}^{OP}(f)\triangleq\int_{0}^{t}f(X_{s},Y_{s})\hbox{d}s$ where $(Z_{s},Y_{s})$ solves the elasto-plastic problem and $(X_{s},Y_{s})$ solves the obstacle problem. The case ($E_{5}$) corresponds to the case $f(z,y)=y$. The case ($E_{5}^{\prime}$) corresponds to the case $f(x,y)=x$. We drop the superscript $(.)^{EP/OP}$ for convenience. In both cases, the process $t\mapsto\Delta_{t}(f)$ has continuous trajectories and admits an asymptotic (in time) variance of the form $\gamma_{f}^{2}t,\>\gamma_{f}^{2}>0$. Thus, with respect to an appropriate scaling in time and space, it is reasonable to expect that it behaves like a Brownian motion in the sense that, taking $(b,T)=(\sqrt{p}\,\beta,p\Theta),\>p$ large, the following approximation can be made

[TABLE]

where

[TABLE]

For an elasto-perfectly-plastic oscillator excited by white noise, this Brownian behavior was already proposed in [13] as a heuristic argument and, recently, a functional central limit theorem (FCLT) has been given in [21]. In the case of a penalization of variational inequalities, a FCLT has been obtained in a general framework including penalization of (3) and (5) (with purely elastic impacts, $e=1$) [33]. The computation of a quantity of type $D$ relies on the computation of quantity of type $C^{\prime}$ (shown in Section 4). See Figure 2 for illustration.

Computational savings

In order to compute the results shown in Figures 2a and b, we used MATLAB codes (running on Intel Xeon E5 2.3Ghz processes) for both our proposed method and the Monte-Carlo method. Our method typically took around one minute, whereas the Monte-Carlo method took around 5 hours to achieve similar accuracy.

5. Numerical computation of $A,B,C$ and $A^{\prime},B^{\prime},C^{\prime}$

In this section, we first explain the probabilistic numerical approach for SVIs of the elasto-plastic (3) and obstacle (5) problems. Then, we present a numerical approach for solving the PDEs with non-standard boundary condition related to the FKf. The method is first presented and applied to the elasto-plastic problem (11), (12) and (13). Then, it is applied to the obstacle problem, (14), (15) and (16).

5.1. Discretization of SVIs and Monte Carlo approach

For the probabilistic numerical scheme of $\{(Z_{t},Y_{t}),\>t\geq 0\}$ of (3), we use a time step $\delta t>0$. Then we construct random variables $\{(Z_{n},Y_{n}),1\leq n\leq N_{\delta t}\}$ and a partition of $[0,T]$, $0=t_{0}<t_{1}<t_{2}<\ldots<t_{N_{\delta t}}=T$. Here, for each $1\leq n\leq N_{\delta t},(Z_{n},Y_{n})$ is an approximation of $(Z_{t_{n}},Y_{t_{n}})$. We refer to [5] for a detailed presentation of the algorithm. To be concise, we only explain the computation of quantities of type $A$. For any choice of well-behaved functions $f$ and $g$, we proceed with the following approximation

[TABLE]

where $\{(Z^{m},Y^{m}),m=1,\ldots,M\}$ is an i.i.d. sequence of trajectories produced by the algorithm.

Details of the probabilistic numerical scheme to compute $\{(X_{t},Y_{t}),\>t\geq 0\}$ in (5) are given in Appendix A. Then we use a similar approximation to that given by (18) in terms of $(X,Y)$ of (5).

5.2. Discretization of PDE problems

To numerically approximate the solutions of (11), (12) and (13), we use a finite difference scheme. We truncate the unbounded domain $D$ to obtain $D_{Y}\triangleq(-1,1)\times(-Y,Y)$, where $Y$ is chosen sufficiently large that the probability of finding the underlying process outside $D_{Y}$ is negligible. We apply a homogeneous Neumann boundary condition at $y=\pm Y$ . We consider a two-dimensional rectangular finite difference grid,

[TABLE]

where $\delta z\triangleq\frac{2}{I-1},\delta y\triangleq\frac{2Y}{J-1}$. Here, $I,J$ are odd integers of the form $2\tilde{I}+1,2\tilde{J}+1$.

The total number of nodes in $\mathcal{G}$ is $N=IJ$. The numerical approximations of $u(z_{i},y_{j},t_{n})$, $v(z_{i},y_{j},t_{n})$ and $w(z_{i},y_{j},t_{n})$ are denoted by ${u}_{i,j}^{n}$, ${v}_{i,j}^{n}$ and ${w}_{i,j}^{n}$ and the corresponding vectors collecting all the unknowns are $\boldsymbol{u}^{n}$, $\boldsymbol{v}^{n}$ and $\boldsymbol{w}^{n}$. We use the notation $f_{i,j},g_{i,j},\psi_{i,j},\varphi_{i,j}$ for $f(x_{i},y_{j}),g(x_{i},y_{j}),\psi(x_{i},y_{j})$ and $\varphi(x_{i},y_{j})$ and the corresponding vectors are $\boldsymbol{f}$, $\boldsymbol{g}$, $\boldsymbol{\varphi}$, $\boldsymbol{\psi}$. Here, $t_{n}\triangleq n\delta t$ discretizes the time and $N_{T}\delta t=T$, $N_{T+h}\delta t=T+h$. Using an implicit Euler method to discretize in time together with finite differences in space, the first conditions in (11), (12) and (13) at the black points in Figure 3 result in

[TABLE]

where $\left(D_{y}\boldsymbol{u}\right)_{i,j}$ and $(L\boldsymbol{u})_{i,j}$ are centered finite differences,

[TABLE]

and first order upwind differences

[TABLE]

with

[TABLE]

The nonstandard boundary conditions (second condition) in (11), (12) and (13) at the red points in Figure 3 are discretized by the same formulae with $L$ replaced by $L_{\pm}$. Here, $(L_{\pm}\boldsymbol{u})_{i,j}$ are defined by

[TABLE]

and

[TABLE]

The Neumann boundary conditions at the points shown in grey results in

[TABLE]

This results in the following linear system to be solved in each time step:

[TABLE]

where $M$ is a sparse $N\times N$-matrix that does not depend on $n$.

For the computational results presented in the remainder of this paper, we use a C code that implements a Monte Carlo (MC) probabilistic simulation to approximate the solution of (3). We use a MATLAB implementation for the PDE approach (23). Here, an LU factorization is used for solving the linear systems arising from the PDE-approach. For time-dependent PDE simulations, the LU factors are computed once upfront and reused throughout the simulation. Implementations are available upon request.

5.3. Numerical results for the elasto-plastic problem

We first present an empirical study on the convergence of the probability of the plastic state ($E_{1}$) and the mean kinetic energy ($E_{2}$). The results presented in Figure 4 provide insight into the dependence of the PDE solution on the domain truncation $L$ and on the number of discretization points. Note in particular, that the convergence of $\mathbb{E}(Y_{T}^{2})$ requires a sufficiently large value of $L$, which is due to the fact that this quantity involves squared $y$-values. Moreover, we also observe that a sufficiently fine mesh with $\tilde{I}=\tilde{J}\geq 100$ is required to properly resolve the PDE problem.

Next, in Figure 5, we present systematic numerical comparisons between the PDE and the probablistic Monte Carlo (MC) approach. In particular, we show comparisons for the quantities ($E_{1}$), ($E_{2}$), ($E_{3}$), ($E_{4}$) and ($E_{5}$). We have chosen discretization parameters, for which we found small approximation errors in our previous tests summarized in Figure 4. As can be seen, the results found from the PDE solution closely track the results from the probabilistic simulation. However, some quantities are harder to approximate than others and it appears that for the correlation structure in the kinetic energy, as shown in Figure 5(d), there is a slight discrepancy between the PDE and probabilistic results. This may be related to the fact that the quantity is quartic in the velocity. From Figure 5, we observe as expected from theory (see [10]) that $(Z_{T},Y_{T})$ has a unique invariant probability measure. Thus, when $T$ becomes large, for any well-behaved function $f$, $\mathbb{E}f(Z_{T},Y_{T})$ becomes constant and $\mathbb{E}\left(\int_{0}^{T}f(Z_{t},Y_{t})\hbox{d}t\right)^{2}$ obeys linear growth with respect to $T$. To the best of our knowledge, except quantities of type $A$ for large time [5, 4], quantities of type $A,B,C$ or $A^{\prime},B^{\prime},C^{\prime}$ (both transient and stationary case) have not been computed using PDE formulations in such a general framework previously.

5.4. Numerical results for the obstacle problem

Next, we present numerical results for the stochastic obstacle problem. As above, for the PDE formulation (14), (15) and (16) we use a finite-difference scheme in space and the implicit Euler method in time. The main difference compared to the discretization for the elasto-plastic problem presented in Section 5.2 is in the discretization of the non-standard boundary conditions to be satisfied at the red points in Figure 3. We recall that the boundary conditions on $D_{T}^{\pm}$ in the PDEs for the elasto-plastic and obstacle problems are different because they reflect different boundary behaviors for the underlying stochastic processes. In the elasto-plastic problem, the boundary condition (plastic phases) solves a boundary value problem whose boundary data is also a part of the problem. In the obstacle problem, the boundary condition (impacts with the obstacle) consists in identifying the values of the solution on $D_{T}^{\pm}$ to those on $x=\pm 1,\mp y>0$.

[TABLE]

where

[TABLE]

We note that $y_{j}$ is defined to be on the grid, but, in general, the quantity $ey_{j}$ will not correspond to a grid point for $e\in(0,1)$. Consequently, the value imposed in $u_{1,j}$ (or in $u_{I,j}$) is a value interpolated between the values of $u$ at the two nearest neighbors of $-ey_{j}$ in the grid. Here we use the notation $[a]$ for the integer part of $a$. In the purely elastic case $e=1$, $j_{1}=J-(j-1)$ and $c_{j}=1$ whereas in the purely inelastic case $e=0$, $j_{0}\equiv\tilde{J}+1$ and $c_{j}=1$.

In Figure 6, we present numerical comparisons between the PDE and probabilistic MC approaches for ($E_{2}^{\prime}$) and ($E_{3}^{\prime}$). As can be seen, the solution of the PDE approach agrees well with the solution of the probabilistic simulation. We also observe that when $T$ becomes large, the expectation (left plot) becomes constant and the variance (right plot) grows linearly. While we expect that ergodicity holds for the obstacle problem, to the best of our knowledge this has not been proven in the literature. Additionally, as far as we are aware of, quantities of type $A,B,C$ or $A^{\prime},B^{\prime},C^{\prime}$ (both transient and stationary case) have not been computed using PDEs in such a general framework before. Observe that ($E_{2}^{\prime}$) and ($E_{3}^{\prime}$) are of type $A$ and $A^{\prime}$ when $T$ is finite.

5.5. Long-time behavior

Next, we present results for the solution of the stationary versions of (11) combined with (13), and of (14) combined with (16) where $\lambda$ is chosen small enough. We remind the reader that $\lambda u$ approximates a quantity of the form $\lim_{t\to\infty}\mathbb{E}f(Z_{t},Y_{t})$ when $\lambda$ is small. We use a similar finite-difference procedure as above for the time-dependent case to obtain a linear system of the form

[TABLE]

Again, we first compare results computed using the discretized PDE approach (24) with results based on Monte Carlo simulations. For the elasto-plastic PDE discretization, we use $\tilde{I}=\tilde{J}=800$, i.e., the overall number of spatial unknowns is $(2\tilde{I}+1)(2\tilde{J}+1)=2,563,201$, and we choose $\lambda=10^{-3}$. For the Monte Carlo simulation, we use $\delta t=10^{-5}$. We present the comparison between the PDE and the Monte Carlo (MC) approach in Table 2. For the obstacle problem, the spatial discretization uses $\tilde{I}=\tilde{J}=500$, amounting to overall $1,002,001$ unknowns. We again use $\lambda=10^{-3}$. Results of the comparison between the PDE and the Monte Carlo approach are given in Table 3. As can be seen, for both the obstacle as well as the elasto-plastic problem, the results obtained with Monte Carlo are close to the PDE results. To study the interplay between the time discretization and the Monte Carlo errors, we considered the case of the growth rate of the variance related to the plastic deformation in the elasto-plastic problem. Approximations of ($E_{1}$) and ($E_{3}$) are shown in Table 1 for $P_{Y}=0.25$ and for different values of $\delta t$ and $T$. The approximation of ($E_{1}$) is given by (18) with $g=0$ and $f(y,z)=\mathbf{1}_{\{|z|=P_{Y}\}}$ and the approximation of ($E_{3}$) is given by

[TABLE]

where $g(y,z)=y\mathbf{1}_{\{|z|=1\}}$. In both approximations, $(Z_{n}^{m},Y_{n}^{m})$ satisfies the probabilistic numerical scheme for (3) and the Monte Carlo sample size is $10^{5}$. For the approximation of ($E_{1}$) we observe a relatively fast convergence towards a constant with respect to $T\geq 3$ whereas for the approximation of ($E_{3}$) a sufficiently large value of $T\geq 150$ is required to see the convergence. For both cases, a sufficiently small value of $\delta t=10^{-4}$ is required. Empirically, the data indicates that the error is roughly halved as the time step is divided by 10.

6. Comparison with [5, 8, 4] for a white noise EPPO at large time

In this section, we compare the approach proposed here with previous techniques employed for a white noise EPPO.

6.1. Approach in [5] for the elasto-plastic problem

If the values of the function $u$ were known at $z=P_{Y},y=0^{+}$ and $z=-P_{Y},y=0^{-}$, the stationary version of (11) would be a degenerate elliptic PDE with a standard Dirichlet boundary condition. This PDE could be solved by a standard numerical approach. However, the challenge in solving this problem resides in the fact that these values are not input data of the problem but part of the solution. This makes it a non-standard problem and more challenging to solve. As a remedy, a superposition approach has been proposed in [5]. To ensure continuity of $u$ at the points $(P_{Y},0)$ and $(-P_{Y},0)$, the linearity of $\eqref{uep}$ allows to compute the solution as a linear combination of the following three local problems:

[TABLE]

with $v(P_{Y},0)=0,v(-P_{Y},0)=0$,

[TABLE]

with $\pi^{+}(P_{Y},0)=1,\pi^{+}(-P_{Y},0)=0$,

[TABLE]

with $\pi^{-}(P_{Y},0)=0,\pi^{-}(-P_{Y},0)=1$.

Then, finding a continuous solution $u$ amounts to finding scalar values $u_{+}$ and $u_{-}$ such that $u=v+u_{+}\pi^{+}+u_{-}\pi^{-}$ is continuous in $(-P_{Y},0)$ and $(P_{Y},0)$, which requires to solve the following linear $2\times 2$ system:

[TABLE]

where

[TABLE]

Remark 2**.**

It is important to note that this technique in [5] is based on the superposition of local PDEs and is thus specific to the case of a SVI modelling an elasto-perfectly-plastic oscillator (EPPO) excited by white noise. It has no natural extension to more general EP or obstacle problems. Indeed, as explained above, solving the KE with non-standard boundary conditions corresponding to an EPPO excited by white noise reduces to solving two linear equations for two unknown scalars. For more general problems as targeted here, one must find unknown functions, and thus the superposition method cannot be employed straightforwardly. This superposition technique has a probabilistic interpretation in terms of novel notion of short cycles as defined in [7].

6.2. Approach in [4, 8] for the growth rate of the variance of $\Delta_{t}$ at large time

In [8], the growth rate of the variance of $\Delta_{t}$ at large time has been characterized. Relying on (3), the authors proposed a novel and simple formulation of the evolution of the system in terms of stopping times in order to identify a repeating pattern in the trajectory, namely long cycles. This concept is summarized next.

Definition 3**.**

A long cycle of the solution $(Z_{t},Y_{t})$ of (3) is a path, enclosed by the stopping times $\tau_{0}$ and $\tau_{1}$ defined below, starting and ending in one of the two points $\{(-P_{Y},0),(P_{Y},0)\}$ which has touched the other point at least once. Similarly, a half long cycle is a path enclosed by the stopping times $\tau_{0}$ and $s_{0}$, or by $s_{0}$ and $\tau_{1}$. In a recursive way, a sequence of stopping times $\tau_{n}$ can be defined where $\tau_{n}$ is the time when the $n$-th long cycle ends.

With the notation

[TABLE]

$\delta\triangleq\textup{sign}(Z_{\tau_{0}})$, which labels the first boundary reached by the process $(Z_{t},Y_{t})$, and

[TABLE]

they obtained a probabilistic expression for the coefficient of the growth rate of the variance of $\Delta_{t}$ as follows:

[TABLE]

where

[TABLE]

In [4], using a Fourier transform approach, an analytic framework for $\mu$ and $\gamma^{2}$ has been proposed.

Remark 4**.**

It is important to note that the technique in [8, 4] of splitting the trajectory in terms of long cycle (an identically independent distributed repeating pattern) is specific to the case of a SVI modelling an EPPO excited by white noise. It has no natural extension to more general EP, obstacle or colored noise problems. Moreover, it cannot be employed for short durations.

7. Conclusions and perspectives

PDEs with non-standard boundary conditions for FKf that describe non-smooth stochastic processes have been (formally) derived and numerically solved. Important examples from engineering have been successfully treated. The present work will pave the way for promising new research. Indeed, the technique presented in this paper is straightforwardly generalizable to a broad range of cases where $F$ and $G$ are time-dependent and to problems driven by colored noise. This framework can also be extended to study Power Spectral Densities by using a system of PDEs with non-standard boundary conditions. Indeed, considering for instance the impact problem, if a complex valued function $u\in C^{\star}(\mathbb{R}^{2}\times[0,T])$ and a real valued function $v\in C^{\star}(\mathbb{R}^{2}\times[0,T])$ satisfy

[TABLE]

together with the boundary conditions for the impact problem, then we have

[TABLE]

Here, $L$ is the infinitesimal generator of $(X,Y)$ away from the boundary.

Appendix A Probabilistic simulation

Below a detailed implementation of the probabilistic simulation for (5). For $T>0,N\in\mathbb{N}$ and $\delta t\triangleq\frac{T}{N}$. We set $\Sigma\in\mathbb{R}^{2\times 2}$ such that

[TABLE]

where

[TABLE]

Here, $\omega\triangleq\dfrac{\sqrt{4k-c_{0}^{2}}}{2}$, it is assumed that $4k>c_{0}^{2}$. For every $(x,y)\in\mathbb{R}^{2}$, we define

[TABLE]

where

[TABLE]

We can now detail the probabilistic simulation for (5), where the main challenge is the incorporation of the obstacles in the simulation. Let $(G_{n,m})_{n=0..N,m=1,2}$ be random independent Gaussian $\mathcal{N}(0,1)$ variables. To compute the $(n+1)$th time step, we attempt to perform the explicit step

[TABLE]

If we find that the $(n+1)$st point does not satisfy the obstacle conditions, for instance since ${X}_{n+1}>1$, we adjust the step length to $\theta_{n+1}\delta t$, with $\theta_{n+1}\triangleq({1-X_{n}})/({{X}_{n+1}-X_{n}})$, and set

[TABLE]

An analogous reduction of the step length is used if the full step computed in (28) satisfies ${X}_{n+1}<-1$.

Appendix B Proofs

B.1. Bridge from solution of (1) to PDEs (when $\mathbf{H}=0$) related to $A,B,C$ and $A^{\prime},B^{\prime},C^{\prime}$.

This is classic but we put it in Appendix for convenience of the reader.

Proposition 5**.**

For any function $\phi\in C^{\star}(\mathbb{R}^{2}\times[0,\infty])$ satisfying $(\mathcal{A}(\lambda,\phi)$,

[TABLE]

the process

[TABLE]

satisfies

[TABLE]

where $L$ is the differential operator

[TABLE]

Thus $M_{t}^{\lambda,\phi}$ is a martingale under $\mathcal{F}_{t}\triangleq\sigma\{W_{s},0\leq s\leq t\}$. Moreover, for any function $\phi^{\prime}\in C^{\star}(\mathbb{R}^{2}\times[0,\infty])$ satisfying $(\mathcal{A}(\mu,\phi^{\prime}))$

[TABLE]

B.2. Proof of Proposition 5

Proof.

From the assumption on $\phi$ the stochastic integral $\int_{0}^{t}e^{-\lambda\tau}\dfrac{\partial\phi}{\partial y}(X_{\tau},Y_{\tau},\tau)\hbox{d}W_{\tau}$ is well defined. Thus, (30) is obtained by using Ito’s formula,

[TABLE]

Therefore $M_{t}^{\lambda,\phi}$ is a martingale with respect to $\mathcal{F}_{t}$. Similarly for $\phi^{\prime}$,

[TABLE]

Hence, (31) is obtained using Ito’s isometry. ∎

B.3. Proof for the backward-in-time parabolic problems.

Proof.

Applying (29)-(30) of Proposition 5 to $\phi=u$, the solution of (7), we get

[TABLE]

Applying (31) of Proposition 5 to $\phi=u$ and $\phi=v$, we get

[TABLE]

Finally, applying (29)-(30) of Proposition 5 to $\phi=w$, the solution of (8), we get

[TABLE]

Therefore,

[TABLE]

∎

B.4. Proof for the degenerate elliptic problems

Proof.

Under Assumption ($\mathcal{A}_{1}(\lambda,\phi)$), using Proposition 5,

[TABLE]

and under Assumption ($\mathcal{A}_{2}(\lambda,\phi)$), direct calculations yields

[TABLE]

∎

Acknowledgement

We thank the editor and the two anonymous reviewers for their helpful comments leading to the improvement of the present manuscript. LM expresses his sincere gratitude to the Courant Institute for being supported as Courant Instructor in 2014 and 2015, when this work was initiated. LM is also supported by a faculty discretionary fund from NYU Shanghai and the National Natural Science Foundation of China, Research Fund for International Young Scientists under the project #1161101053 entitled “Computational methods for non-smooth dynamical systems excited by random forces” and the Young Scientist Program under the project #11601335 entitled “Stochastic Control Method in Probabilistic Engineering Mechanics”. LM also thanks the Department of Mathematics of the City University of Hong Kong for the hospitality. JW acknowledges support from the SAR Hong Kong grant [CityU 11306115] “Dynamics of Noise-Driven Inelastic Particle Systems”.