1. No category

Proceedings of the 9th International Conference on Structural Dynamics, EURODYN 2014
Porto, Portugal, 30 June - 2 July 2014
A. Cunha, E. Caetano, P. Ribeiro, G. Müller (eds.)
ISSN: 2311-9020; ISBN: 978-972-752-165-4
Numerical analysis of stochastic differential equations
arising in dynamics of structures
1
Pierre Bernard1 , Frédéric Bernardin2 , Jun Shao3
Laboratoire de Mathématiques, UMR CNRS/UBP 6620, Blaise Pascal University, Campus des Cézeaux, 63171,
Aubière CEDEX, France
2
CEREMA, Direction territoriale Centre-Est, Département Laboratoire de Clermont-Ferrand, 63017,
Clermont-Ferrand CEDEX, France
3
Institut Pascal, UMR CNRS/UBP 6602, Campus des Cézeaux, 63171, Aubière CEDEX, France
email : [email protected], [email protected], [email protected]
ABSTRACT : There is a huge litterature dedicated to the numerical analysis of stochastic differential equations (SDE).
The SDEs encountered in structural dynamics have some particularities. The drift term is very important and needs
specific schemes. Structures are in most cases slightly dissipative systems. It is important that the schemes used do not
present the following drawbacks : either destroy the mechanical energy dissipation by introducing a numerical energy
increase, or introduce some numerical energy dissipation of the same order (or higher) than the physical dissipation.
Moreover, the excitation can be modelled by white noïse or coloured noïse. In the first case, we deal with stochastic
differential equations in Ito’s or Stratonovich sense. In the second case, the excitation has more regular trajectories and
in many cases we can use schemes from the analysis of ordinary differential equations. That is why we will bring much
attention to schemes with conservative properties when applied to physical conservative systems. This will give rise
to what we shall call symplectic numerical schemes for SDEs. We will also be concerned by numerical algorithms for
multivalued SDEs, or stochastic variational inequations, appearing for example in the analysis of dynamical response
of nonlinear structures under seismic excitation (elasto-plastic or hysteretic behavior under random loads).
KEY WORDS : numerical schemes ; ordinary differential equations ; stochastic differential equations ; symplectic
schemes.
1
Some schemes for ODE’s issued from mechanics
Ordinary Differential Equations (ODEs) issued from
dynamics are considered in this section, that is to say second order equations equations of the form
Newmark method
The following relations are used :
∆x˙ i
∆xi
= h¨
xi + γh∆¨
xi
2
h
¨i + βh2 ∆¨
xi
= hx˙ i + x
2
(2)
x
¨(t) = f (x, x(t))
˙
+ p(t)
(1) where ∆x = x
i
i+1 − xi for a sequence xi , and 0 < γ < 1 et
with initial conditions x(0) = x0 , x(0)
˙
= x˙ 0 . Some regu- 0 < β ≤ 1/2γ are parameters. The most popular choice is
γ = 1/2 and 1/6 ≤ β ≤ 1/4.
larity is assumed for f and p.
For nonlinear oscillators, it willl be necessary to comThere are two possibilities : either keep second order
derivatives, either change to phase space, with position- pute a tangent stiffness at each step, or to use an iterative
velocity variables, and associate a dynamical system with procedure (see [4]).
firts order derivatives. That is go from Lagrange equations
to Hamilton formulation. When cinetic energy is a quadra- 1.2 Schemes for dynamical systems
tic form in the velocity, the two approaches are equivalent.
First order ODE’s are considered in this paragraph.
As it will be necessary to solve equations with many
trajectories of stochastic processe as inputs, we will concen- Note that one can allways switch in this situation by extrate our attention on schemes which do not need too large pessing the dynamics in the phase space. The dynamical
system defined by the ODE
computations.
1.1
Schemes with second order derivatives
y(t)
˙ = f (y(t)),
(3)
where y ∈ IR2n , defines the continuous flow φt , t ∈ [0, T ]
caracterised by φt (y0 ) = y(t) for all t ∈ [0, T ], y0 ∈ IR2n .
For the central difference method with step h , x
¨i
To the numerical scheme resolving (3) is associated the
is approximated by (xi+1 − 2xi + xi−1 )/h2 and x˙ i by discrete flow Φnh such that Φh (y0 ) = y¯1 , where y¯1 is the nu(xi+1 −xi−1 )/2h. These approximations are of second order merical solution at step h issued from y0 computed using
in h. Problems of stability can arise. The stability analysis the scheme. The ODE (3) can be rewriten under the inteis not simple, we cannot treat this question here. One just gral form
Z t
has to know that he must be very carefull with the choice
y(t)
=
y
+
f (y(s))ds
(4)
0
of the time step.
Central difference method
O
2781
Proceedings of the 9th International Conference on Structural Dynamics, EURODYN 2014
and the first Euler schemes are obtained using elementary 1.3 Symplectic schemes
approximations of the integral :
Moreover, as we deal with mechanical systems, an important property of the flow associated to Hamilton EquaExplicit Euler scheme
tions is symplecticity, that is to say the conservation of the
The integral on [0, h] is approximated by hf (y0 ). So
Hamiltonian form of the system. Using a numerical scheme
to solve an Hamilton Equation, it is important that the disy1 = y0 + hf (y0 )
(5)
crete flow Φh be symplectic.
This scheme is explicit. knowing, y0 , y1 is immediately calculated.
Symplectic flows
Implicit Euler scheme
The interested reader can refer to [1]. Concerning the
domain
of numerical methods, the reference here is [7].
The integral on [0, h] is approximated by hf (y1 ). So
Let J denote the 2d-dimensional symplectic matrix, and
y1 = y0 + hf (y1 )
(6) ω the symplectic form of the same dimension :
This scheme is implicit : one has to solve an implicit equa"
#
0
I
tion to compute y1 .
ω(ξ, η) = ξ T Jη; where J =
(10)
−I 0
Mean point Euler scheme
The integral on [0, h] is approximated by hf ((y0 +
y1 )/2). So
Definition 1 A linear application A : IR2d → IR2d will
y1 = y0 + hf ((y0 + y1 )/2)
(7) be called symplectic if it leaves invariant the symplectic
form : ω(Aξ, Aη) = ω(ξ, η) ∀ξ, η ∈ IR2d . or equivalently
This is also implicit.
if AT JA = J. Such an application is also called canonical
transform.
Trapezoidal Euler scheme
A differentiable mapping g : U → IR2d where U is an
The integral on [0, h] is approximated by h(f ((y0 ) + open set in IR2d is symplectic if its differential is.
f (y1 ))/2. So
y1 = y0 + h(f ((y0 ) + f (y1 ))/2
(8) This property can be interpreted as the conservation of
bidimensional oriented areas.
This is also implicit. These schemes are not very different
Let us now consider an Hamiltonian system in IR2d defifrom one another. However, they have different stability ned by the Hamiltonian H(p, q) using usual notations. The
properties.
associated dynamics is given by the system of Hamilton
equations :
Runge Kutta methods
y˙ = J −1 ∇H(y).
(11)
The principal is to use evaluations of f (y(t)) at some inThe following result was proved by Poincaré :
termediate points.The most popular among Runge-Kutta
(RK) methods is the RK4, which is as follows :
Theorem 1 If H(p, q) is a C 2 -Hamiltonian on the open
h
set U ⊂ IR2d and φt the associated flow, then φt is sym(9)
yn+1 = yn + (K1 + 2k2 + 2k3 + k4 )
plectic.
6
h
k1 = f (yn ) ; k2 = f (yn + k1 ) ;
The advantage of symplectic applications is to trans2
form
Hamiltonian systems into Hamiltonian systems.
h
k3 = f (yn + k2 ) ; k3 = f (yn + hk3 )
2
2 Let ψ : U → V a coordinate change that
This is an explicit method. Stability is conditional. It is a Theorem
1
is
a
C
-diffeomorphism
Il ψ is symplectic, any Hamiltofourth order method, under conditions of regularity for f .
nian system y˙ = J −1 ∇H(y) becomes in the new variables
z = ψ(y) the Hamiltonian system z˙ = J −1 ∇K(z) where
Collocation methods
K(z) = H(y). Conversely, if ψ transforms any HamiltoThe idea is to choose a finite dimensional space of po- nian system in Hamiltonian system, ψ is symplectic.
lynomials candidate solutions and a number of points, the
Let us now consider the associated numerical schemes.
collocation points, and to select that polynomial which is
an exact solution at the collocation points. It is a special class of implicit Runge Kutta methods. The trapezoidal Euler method is a particuler case of collocation me- Definition 2 A one step numerical method is called symthod. The Gauss-Legendre method use the points of Gauss- plectic if the application y1 = Φh (y0 ) is symplectic when
applied to a regular Hamiltonian system.
Legendre quadrature as collocation points.
2782
Proceedings of the 9th International Conference on Structural Dynamics, EURODYN 2014
2
Some symplectic schemes
Remark 1 Explicit, implicit and trapezoidal Euler
schemes are not symplectic ; the mean point Euler method
scheme is symplectic.
Partitioned systems
Main schemes for SDE’s issued from Mechanics
The main references here are the book [8] and the review article [9]. The articles [3] and [6] contains some analysis of schemes more specific for Mechanics.
Let us assume that the variables are partitioned (for
2.1
example in the phase space : position and velocity) so that
the differential system can be written
u˙ = a(u, v)
v˙ = b(u, v)
(12)
In that case the partitioned Euler schemes can be considered to solve such a system :
un+1 = un + ha(un , vn+1 )
vn+1 = vn + hb(n , vn+1 )
(13)
un+1 = un + ha(un+1 , vn )
vn+1 = vn + hb(un+1 , vn )
(14)
or
where one variable is treated by the implicit and the other
by the explicit Euler method. This method is symplectic.
Some basic definitions and results
concerning the stochastic case
We will have to distinguish the EDS in the Itô’s sense
(EDSI) and in the Stratonovich sense (EDSS). The second
type is more natural to modelize mechanical problems with
random excitation of White Noise type, but the mathematical analysis uses Itô’s acception. The definition of weak
and strong solutions, then weak and strong convergence of
numerical schemes, and of order of convergence, are defined in the given reference. As well as the trajectorial unicity and unicity in law of solutions. To give a rough idea,
strong refers to trajectories, while weak refer to probability laws. These last are generally the interesting objects
for mechanics since the goal is to compute probabilities.
We will denote by d the differential in the Itô sense, and
by δ the differentials in Stratonovich sense. Let us recall
that the equation :
dX(t) = b(X(t))dt + σ(X(t))dW (t); X(0) = X0
(16)
One exemple : the Lotka-Volterra problem
where W is a Brownian motion is a way of writing the
This model is the prey-predator dynamics issued from integral equation
Z t
Z t
biology. But it is a good test for algorithms. It models the
evolution in time of the number of predators u(t) and the
X(t) = X0 +
b(X(s))ds +
σ(X(s))dW (s) (17)
0
0
number of preys v(t) at time t.
u˙ = u(v − a)
v˙ = v(b − u)
(15)
in which the first integral is a classical lebesque integral
while the second one is an Itô integral. It is the same for
EDS in the sense of Stratonovich, replacing dW by δW .
The numerical schemes will rest on the approximation
of the integrals and on stochastic Taylor formulas obtained by iteration of Itô’s formula. We need to approximate
first the deterministic integral, then the stochastic integral.
The question of implicit schemes will be perturbated by the
non-anticipative terms in stochastic calculus.
where a and b are positive constants.
This systems admits an invariant (Lagrangian)
Φ(u, v) = bLn(u) − u + aLn(v) − v. Let us solve this system
using explicit Euler, implicit Euler, mean point Euler and
Trapeze Euler.
The flows associated to these different schemes are represented in figure 1, with time step 0.1 and duration
Euler Maruyama scheme
T = 10.
The time interval is discretized at step h on the interval
[0, T ]. This scheme is like the Euler scheme in the deterministic case, and will be called Euler scheme in the rest of
the text
6
Explicit Euler
Implicit Euler
Mean point
Trapezoidal
5
4
3
Y0 = X0 , Yn+1 = Yn + hb(Yn ) + σ(Yn )∆Wn
(18)
2
where ∆Wn = W ((n + 1)h) − W (nh) is the sequence of
accretions of the Brownian motion, that is a sequence of
independant zero mean Gaussian variables with variance
h. Hence, for the simulation, ∆Wn will be substituted by
√
Figure 1. Flows associated to Lotka-Voltera problem obhNn where (Nn ) is a sequence of reduced independant
tained by different schemes
Gaussian r.v.
1
0
0
0.5
1
1.5
2
2.5
3
3.5
4
4.5
5
2783
Proceedings of the 9th International Conference on Structural Dynamics, EURODYN 2014
Continuous approximation
Order of strong convergence of the Euler scheme
We defined an approximation of a continuous process
Let us assume that the functions b and σ are C 2 with
by a discrete one. An approximation by another continuous bounded derivatives of first and second order. Then, the
process will sometimes be needed. This can be obtained by Euler scheme is of strong order 12 .
linear interpolation on each interval :
Y (t) = Y + (t/h − n)(Y
− Y ), t ∈ [nh, (n + 1)h] (19) Order of weak convergence of the Euler scheme
n
n+1
n
Under some conditions (see [8]) on the coefficients, the
Another method will be obtained by freezing in the terms
Euler
scheme is of weak order 1.
of (16) the coefficients on each interval at the value at the
Remark
: to obtain a weakly convergent scheme, it is
begining of the interval. Denote btcn = Ent(nt)/n where
not
necessary
to use Gaussian variables ∆Wn . One can use
Ent is the integer part. The approximation then is solution
a
random
walk
for example.
of the following EDSI :
An order 1/2 can be unsatisfactory. That is why was
Y h (t) = Ynh + b(Ynh )t + σ(Ynh )dW (t), t ∈ [btcn , btcn + h] introduced the following :
(20)
Milshtein scheme
Weak and strong convergence
Firts consider the dimension 1. Using a Taylor developFor the strong approach, we compare the exact solution ment of σ(·) in (16), one obtains for small t :
with its approximation using the same trajectory of the
X(t) ≈ X0 + b(X0 )t + σ(X0 )(W (t) − W (0) (21)
Brownian motion W . This will be a trajectorial approximaZ t
tion. The error will be measured by E(|X(T ) − Y (T )|p )1/p ,
0
+ σ(X0 )σ (X0 )
(W (s) − W (0))dW (s)
generally with p = 1 ou p = 2. One can also use an uniform
0
error on all time steps, like supk=0,...,N E(|X(kh) − Yk |) or
supt∈[0,T ] E(|X(t) − Y (t)|).
The presence of an iterated Itô’s integral makes things comOne can also retain almost sure uniform convergence plicated, since this should be simulated. If W is one dimenon [O, T ].
sional, it is simple since, by Itô’s formula,
For the weak approach, we will compare the laws, for
Z t
1
example using moments. So the error willl be mesasured by
(W (s) − W (0))dW (s) = (W (t)2 − t)
2
|E(X(T )p ) − E(Y (T )p )| for convenient values of p. More
0
generally, one can use |E(g(X(T )) − E(g(Y (T ))|) for some
class of test functions G. If G contains the polynomials, this This gives the Milshtein’s scheme :
will imply a control of the error on moments.
Y(k+1) = Yk + b(Yk )h + σ(Yk )∆k W
1
Order of convergence
+
σ(Yk )σ 0 (Yk )((∆k W )2 − h).
2
In a stochastic context, the notions of convergence are
But in dimension more than one, the extension of the
more involved than in a deterministic context. Let us deformula
(21) is, component by component :
h
note Y the approximation sequence of X(t) resulting from
a numerical scheme with step h. See [8]
X i (t) ≈ X0i + bi (X0 )t + Σdj=1 σji (X0 )(W j (t) − W j (0) (22)
Z t
Definition 3 The scheme Y h is strongly convergent
+ Σdj,l=1 ∂σji (X0 )σli (X0 )
(W l (s) − W l (0))dW j (s)
with order γ > 0 if there exists a constant C > 0, which
0
does not depend on h and h0 > 0 such that
The multiple Itô’s integrals present here introduce some
Err2 (h) = [E(|X(T ) − Y h (T )|2 )]1/2 ≤ Chγ
difficulties, since they are non continuous fonctionals of the
for all h ∈ [0, h0 [.
Brownian motion W , which will forbid trajectorial approxiThe scheme Y h is weakly convergent with order mations, and moreover there law is difficult or impossible
β relatively to the class of functions G if, for any g ∈ G, to express with Gaussian distributions.
there exists a constant C(g, T ) depending on g but not on
h, and h0 > 0 such that
Order of strong convergence of the Milshtein’s
err1,g,T (h) = |E(g(X(T )) − E(g(Y h (T )))| ≤ C(g, T )hβ
scheme
Assume b and σ are C 2 and the two firts order derivatives bounded. Then the Milshtein’s scheme is of strong
One can also consider an order of almost sure conver- order 1, that is there existe a constant C(T ) such that
gence α if there exists an almost surely finite r.v. C and
h0 > 0 déterministic such that
E(|X(T ) − Y h (T )|2 )1/2 ≤ C(T )h.
for all h ∈ [0, h0 [.
supt∈[0,T ] |X(t) − Y h (t)| ≤ Chα for all h ∈ [0, h0 [.
2784
The result is also available with L1 instead of L2 -norms.
Proceedings of the 9th International Conference on Structural Dynamics, EURODYN 2014
Cameron-Clark commutativity condition
Some schemes satisfying these assumptions
– Explicit and implicit Euler scheme
However, there is a situation in which the simulation is
– Milshtein scheme in dimension 1
possible. Let X live in IRd and W in IRm , so that σ(x) is a
– Stochastic Newmark scheme
matrix with d lines and m columns, with coefficients σji (x).
Let us introduce the first order differential operators assoj
Theorem 3 Let T > 0 given. Let X(t) a solution of (16)
ciated to the columns of σ : Lj = Σm
l=1 σl ∂l , j = 1, · · · , d.
and Xn (t) a continuous approximation satifying (26) :
Consider the following commutativity condition :
∀j1 , j2 ∈ {1, · · · , d}, ∀k ∈ {1, · · · , m}, Lj1 σkj2 = Lj2 σkj1
dXn (t)
=
bn (Xn (btcn )) dt + σn (Xn (btcn )) dW (t)
(23)
+
ϕn (Xn (btcn )) dZn (t),
Note that this condition is satisfied in the following cases :
X
(0)
=
x
.
n
n,0
– m = 1,
– σ is constant, independant of x,
Then under some technical asumptions, for any α ∈ 0; 21 ,
– d = m and σ is diagonal,
there exists a random variable η such that, almost surely,
– σ(x) is linear.
sup |X − Xn |(t) ≤ η.n−α .
Proposition 1 Let us assume the condition (23) is satis0≤t≤T
fied. Then the Milshtein scheme can be written :
l
Y(k+1)
l
j
= Ykl + bl (Yk )h + Σm
(24)
j=1 σj (Yk )∆k W
1
Σj <j (Yk )Lj1 σjl 2 (Yk )((∆k W j1 )(∆k W j2 )
+
2 1 2
1 m
+
Σ (Yk )Lj σjl (Yk )((∆k W j )2 − h)
2 j=1
The Milshtein’s scheme is optimal in the following sense
[5] :
Stochastic Newmark scheme
As the Newmark scheme is very popular in engineering,
it is implemented in many codes. However the quality of
this scheme depends on the regularity on the sollicitaion,
that is supposed at least C 3 in classisal studies. Then this
scheme is used with white noïse. It appears then that is
is very smoothing and filters high frequencies too much.
Proposition 2 Under the condition (23), the Milshtein’s See [3]. So in this reference modifications of the Newmark
scheme has the best order of convergence in quadratic mean scheme are presented to use it with white noïse.
among all schemes with fixed time step h using only values
of the Brownian motion W at times kh.
2.2 Stochastic symplectic schemes
Almost sure trajectorial convergence
We refer to the paper by G.Fleury (6) for more results
on this topic. The extension is as follows : the classical assumption that the coefficients should be globally Lipschitz
is not satisfied in many examples (Duffing or Van der Pol
oscillator, and so on). The very general following schemes
are considered :
Xn,k+1
Xn,0
=
Xn,k + bn (Xn,k ).h + σn (Xn,k ).∆n,k W
+
ϕn (Xn,k ).∆n,k Zn ,
= xn,0 ,
These schemes are dedicated to SDE describing the dynamics of systems that can be considered as perturbations
of Hamiltonian systems. The perturbation is double : a dissipative term "not too big" and the action "not too strong"
of the noise. To be more precise on what we mean by "not
too big", see the reference [2]. As we know from [8], use
implicit schemes for the stochastic part leads to important
difficulties.
Mean point scheme
The dynamics is described by the Itô SDE :
dX(t) = b(X(t))dt + σ(X(t))dW (t); X(0) = X0 ; (26)
where the discrete time process (Xn,k )k∈IN is an observawhere
tion at times nk (k ∈ IN) of the solution of :
dXn (t)
= bn (Xn (btcn )) dt + σn (Xn (btcn )) dW (t)
(25)
+ ϕn (Xn (btcn )) dZn (t),
Xn (0) = xn,0 .
where


Zn



 b
n

σn




ϕn
are continuous martingales,
: IRh −→ IRh arevector − valuedfunctions ,
: IRh −→ Mh0 ,h (IR) et
: IRh −→ Mh”,h (IR) are matrix − valued functions.
Moreover, (σn ), (bn ) are supposed to converge to σ and b
respectively, and (ϕn ) are uniformly bounded.
b(x, v) = J −1 ∇H(x, v) + εc(x, v)
!
σ1 (x)
X = (x, v), σ(x, v) =
0
(27)
We propose to use the mean point Euler scheme for the
deterministic part, which is known to be symplectic, and
the Milshtein sheme for the random part. In dimension 1,
the following is obtained :
Yk+1
= Yk + b((Yk + Yk+1 )/2)h + σ(Yk )∆k W (28)
1
+
σ(Yk )σ 0 (Yk )((∆k W )2 − h); Y0 = X0
2
2785
Proceedings of the 9th International Conference on Structural Dynamics, EURODYN 2014
Problems that can arise using implicit schemes for is univalued. As an example, the Coulomb’s friction law
is defined thanks the multivalued operator ACoulomb on R :
the stochastic part

Consider the one dimensional linear Itô’s SDE :
−1 if x < 0


(32)
ACoulomb (x) = [−1, 1] if x = 0

dX(t) = aX(t)dt + bX(t)dW (t); X(0) = X0
(29)

+1 if x > 0
and take the implicit Euler scheme with step h for both
3.2
the deterministic and the stochastic part
Yk+1 = Yk + ahYk+1 + bYk+1 ∆k W
(30)
In this example, one can compute explicitly
Some definitions
MSDE’s
and
results
on
We consider the following MSDE :
(
dXt + A(Xt )dt 3 b(t, Xt )dt + σ(t, Xt )dWt , 0 ≤ t ≤ T,
X0 = x0
Yn =
− ah − b∆k W ).
(33)
Then E[|Yn |] = +∞√if (∆k W ) is a qequence of independant where A is a multivalued and maximal monotone operator
r.v. with law N (0, h).
on Rd . In order to make sense to the inclusion (33), we
√But, using a sequence of independant r.v. taking values introduce a stochastic process (Kt )0≤t≤T a.s. with boun± h with probability 1/2, the r.v. Yn are bounded if h is ded variation on [0, T ], vanishing at 0, satisfaying dKt ∈
A(Xt )dt and such that :
small enough.
X0 Πn−1
k=0 1/(1
dXt = b(t, Xt )dt + σ(t, Xt )dWt − dKt .
3
Multivalued SDE’s
(34)
We give here the result of existence and uniqueness of
a solution to (33) ([10]).
We refer here to [10] for the mathematical statements of
multivalued SDE’s, to the book [11] for the application of Proposition 3 Let us assume that the interior of the doMSDE’s to dynamics of structures and to the book [12] for main of A is a non-empty set :
mechanical applications and numerical aspects of MSDE’s.
We mention here recent works [15] using the framework of
Int(D(A)) 6= ∅;
(35)
stochastic variational inequalites.
Let us furthemore assume that there exists two constants
C(T ) > 0 and γ > 0, such that for all (s, t) ∈ [0, T ]2 and
3.1 Some definitions and results on multi- x, y ∈ Rd ,
valued and maximal monotone operakb(t, x)−b(s, y)k+kσ(t, x)−σ(s, y)k ≤ C(T )(|t−s|γ +kx−yk)
tors
(36)
MSDE’s are used in random mechanics for the modeling of dry friction, elastoplasticity or impacts under random sollicitations. Their drift coefficient are defined thanks
to multivalued operator, that is maps from Rn into P(Rd )
(n ≥ 1) where P(Rd ) denotes the set of parts of Rd ). We
consider here only maximal monotone multivalued operator (see [14]), that is multivalued monotone operators without “hole”. In this subsection we give well-known results
about maximal monotone operator (see [14] for more explanations). The domain D(A) of A : Rd → P(Rd ) is the
set of all x ∈ Rd such that A(x) is not empty. We denote by Gr(A) the graph of A, i.e. Gr(A) = {(x, y) ∈
Rd , x ∈ D(A), y ∈ A(x)}. We say that A is monotone if
for all x, x0 ∈ D(A), and for all y ∈ A(x) and y 0 ∈ A(x0 ),
hx − x0 , y − y 0 i ≥ 0 is valid. For A a multivalued operator,
we define its inverse A−1 or, what is equivalent, Gr A−1 ,
as the graph obtained from Gr A by exchanging the order
of the elements in any pair (x, y) ∈ Gr A. We recall the
following result which will be usefull to define a numerical
scheme for MSDE’s : for λ > 0, the a priori multlivalued
map
Jλ = (I + λA)−1
(31)
2786
et
kb(t, x)k + kσ(t, x)k ≤ C(T )(1 + kxk).
(37)
Then there exists a unique strong solution to (33).
3.3
A generic
MSDE’s
numerical
scheme
for
For the general MSDE (33), an implicit-explicit Euler
scheme can be well-defined. The scheme is defined as follows :
X 0 = x0 ,
∀n ∈ {0, . . . , N − 1},
(38a)
X n+1 − X n + hA(X n+1 ) 3 hb(tn , X n ) + σ(tn , X n ) Wtn+1 − Wtn .
(38b)
By virtue of (31), (38b) is equivalent to :
∀n ∈ {0, . . . , N − 1},
X n+1 = Jh X n + hb(tn , X n ) + σ(tn , X n ) Wtn+1 − Wtn .
(39)
Proceedings of the 9th International Conference on Structural Dynamics, EURODYN 2014
Proposition 4 Let us assume assumptions of Proposition We have :
3. We assume moreover that there exists an integer p ≥ 4
xn+1 − xn
and a constant C > 0 such that :
x˙ n+1 − x˙ n
Ekx0 kp < ∞
(40)
vn+1
and
= hx˙ n+1
(45)
= −chx˙ n+1 − khvn+1 + Σ∆n W
(46)
= max[min(vn + xn+1 − xn , Y ), −Y ](47)
and then afer calculations :
kprojA(x) 0k ≤ C(1 + kxkp/2 ),
∀x ∈ D(A) :
(41)
xn+1
vn+1
= gn (vn+1 )
(48)
= max[min(vn + gn (vn+1 ) − xn , Y ), −Y ](49)
where projA(x) denotes the projection on the convex set
A(x). Then the scheme (38) strongly converges to the so- where
lution of (33) :
gn (x) = 1/(1+ch)[xn (2+ch)−xn−1 −kh2 x+Σ∆n W ] (50)
lim E sup kXi − Xti k2 = 0.
h→0
0≤i≤N
(42) We have to solve the system (48)-(49). Let us denote by x?
the solution of x? = vn + gn (x? ) − xn :
If we assume morever that σ is bounded then :
E sup kXi − Xti k2 = O((h log(1/h))1/2 ),
0≤i≤N
x? = k1 [(2 + ch − kh2 )xn − xn−1 − kh2 vn + Σ∆n W ]
(43) with
k1 = 1/(1 + ch + kh2 )
Proof 1 We refer to [13] for the proof.
3.4
Application to an elastoplastic system
under white Noise
The system is represented by the triplet X(t) =
(x(t), x(t),
˙
v(t)) and is governed by the MSDE in R3 :
dX(t) + A(X(t))dt 3 BX(t)dt + ΣdW (t) where Σ > 0,
W is the real standard Brownian motion and :




0
0 1
0




 , β = σ −1 ,


B=
0

0 −c −k  X, A(X) = 
0
1
0
β(X3 /Y )
If |x? | < Y , then vn+1 = x? and xn+1 = gn (x? ). If x? > Y ,
vn+1 = Y and xn+1 = gn (Y ). Indeed, the map gn is decreasing and then Y < vn + gn (x? ) − xn < vn + gn (Y ) − xn .
It follows that max[min(vn + gn (Y ) − xn , Y ), −Y ] = Y .
If x? < −Y by similar arguments, vn+1 = −Y and
xn+1 = gn (−Y ). The solution of the system (48)-(49) is
given by xn+1 = gn (vn+1 ) and :
 ?
?

 x if |x | < Y
?
Y if x > Y
vn+1 =


−Y if x? < −Y
Algorithm for the scheme with middle point
scheme applied to the linear part
where c, k and Y are positive real constants. We will illus- We discretize the MSDE (38) as follows :
trate in the sequel the scheme defined in (39) and two
schemes for wich the linear part of the drift is discreti- Xn+1 − Xn + hA(Xn+1 )h 3 1 hC(Xn+1 + Xn ) + Σ∆n W
2
zed thanks to an implicit Euler scheme and a mean point
scheme.
We have :
1
xn+1 − xn =
h(x˙ n + x˙ n+1 )
(51)
Algorithm for the scheme (39)
2
1
After calculations, we obtain the following algorithm :
x˙ n+1 − x˙ n = − ch(x˙ n + x˙ n+1 )
2
1
x1 = x0 + hx˙ 0 ,
− kh(vn + vn+1 ) + Σ∆n W
2
v1 = max[min((v0 + x1 − x0 , Y ), −Y ]
vn+1 = max[min(vn + xn+1 − xn , Y ), −Y ]
2
xn+1 = (2 − ch)xn − (1 − ch)xn−1 ) − kh vn−1
The second equation gives :
+Σ∆n−1 W
1
vn+1 = max[min(vn + xn+1 − xn , Y ), −Y ]
x˙ n+1 (1+ch/2)− x˙ n (1−ch/2) = − kh(vn +vn+1 )+Σ∆n W
2
since (I + hβ)−1 (z) = max[min(z, Y ), −Y ].
This last equation and (51) lead to :
h
1 − ch/2
1
Algorithm for the scheme with implicit Euler x˙
(xn+1 − xn ) + Σ∆n W
n+1 = −k (vn + vn+1 ) +
4
h
2
scheme applied to the linear part
This last relation expressed with n instead of n + 1 gives :
We discretize the MSDE (38) as follows :
Xn+1 − Xn + A(Xn+1 )h 3 CXn+1 h + Σ∆n W
h
1 − ch/2
1
(xn − xn−1 ) + Σ∆n−1 W
(44) x˙ n = −k 4 (vn + vn−1 ) +
h
2
2787
Proceedings of the 9th International Conference on Structural Dynamics, EURODYN 2014
References
and, by (51) :
xn+1 − xn
1
h(x˙ n + x˙ n+1 )
2
h2
= −k (vn+1 + 2vn + vn−1 )
8
1 − ch/2
+
(xn+1 − xn−1 )
2
h
+ Σ(∆n−1 W + ∆n W )
4
=
Finally :
xn+1
vn+1
= gn (vn+1 )
(52)
= max[min(vn + gn (vn+1 ) − xn , Y ), −Y ](53)
where
gn (x) =(2/(2 + ch))[2xn − (1 − ch/2)xn−1
h2
(x + 2vn + vn−1 )
(54)
4
h
+ Σ(∆n−1 W + ∆n W )]
2
Denoting by x? the solution of x? = vn + gn (x? ) − xn :
−k
x? = k2 [(2 − ch)(xn − xn−1 ) + (2 + 2ch − kh2 )vn
− (kh2 /2)vn−1 + hΣ(∆n−1 W + ∆n W )]
2
with k2 = 1/(2 + ch + k h2 ).
As in the previous case, the solution of the system (52)(53) is given by xn+1 = gn (vn+1 ) and :
 ?
?

 x if |x | < Y
Y if x? > Y
vn+1 =


−Y if x? < −Y
We illustrate on figure 2 the three schemes introduced
in this section for a time step h = 10−4 with T = 104 ,
c = 0.12, k = 900, Y = 0.04, Σ = 2, x0 = 0.03, v0 = x0 ,
x˙ 0 = 0.
1. V.I.Arnold Mathematical Methods of Classical Mechanics Springer Verlag, New York, second edition,
1989.
2. P.Bernard Stochastic averaging : some methods and
applications IUTAM Symposium on Nonlinear Stochastic Dynamics, p. 29-41. Kluwer. 2003.
3. P.Bernard, G.Fleury, Stochastic Newmark scheme,
Probabilistic Engineering Mechanics 17, 45-61, 2002.
4. A.K.Chopra, Dynamic of Structures, Printice Hall,
1995.
5. J.M.Clark, R.J.Cameron The maximum rate of
convergence of discrete approximations for stochastic differential equations Stochastic Differential Systems Filtering and Control, vol 25 of Lecture Notes
in Control and Information Sciences, Springer, 1980.
6. G. Fleury, Convergence of Schemes for Stochastic
Differential Equations, Probabilistic Engineering Mechanics, 21, 35-43, 2006.
7. E.Hairer, C.Lubich and G.Wanner, Geometric Numerical Integration, Springer, Second edition, 2006.
8. P.E.Kloeden, E.Platen, Numerical Solution of Stochastic Differential Equations, Springer, Third Printing, 1999.
9. D. Talay, Simulation of Stochastic Differential Systems, pp. 54-96. Probabilistic Methods in Applied
Physics, P. Kree & W. Wedig eds, Springer, Lecture
Notes in Physics 451, 1995.
10. E. Cépa, Equations différentielles stochastiques multivoques, pp. 86-107., Lecture Notes in Mathematics,
Séminaire de Probabilités, XXIX, 1995.
0.1
Explicit Euler
Implicit Euler
Mean point
0.05
11. P. Krée and C. Soize, Mathematics of random phenomena : random vibrations of mechanical structures,
D. Reidel Pub. Co., 1986.
0
x
12. J. Bastien, F. Bernardin and C.-H. Lamarque, Non
Smooth Deterministic or Stochastic Discrete Dynamical Systems : Applications to Models with Friction
or Impact, Wiley, 2013.
−0.05
−0.1
−0.15
0
1000
2000
3000
4000
5000
6000
7000
8000
9000
10000
t(s)
Figure 2. An example of a simulated trajectory of x by the
three schemes.
ACKNOWLEDGMENTS
The authors thank the French National Research
Agency for having supported this work.
2788
13. F. Bernardin, Multivaled stochastic differential equations : convergence of a numerical scheme, pp. 393415, Set-Valued Analysis 11, 2003.
14. H. Brézis, Opérateurs maximaux monotones et semigroupes de contractions dans les espaces de Hilbert,
North-Holland Publishing Co., 1973.
15. L. Mertz, Inéquations variationnelles stochastiques et
applications aux vibrations de structures mécaniques,
PhD Thesis, Paris, 2011.