0895-7177188$3.00+ 0.00 Pergamon Press plc
Math1Compur.Modellin~,Vol. 10.No. 3. PP. 207-216.1988 Printed in Great Britain
A STOCHASTIC MODEL IN CONTINUUM MECHANICS: TIME EVOLUTION OF THE PROBABILITY DENSITY IN THE RANDOM INITIAL BOUNDARY-VALUE PROBLEM I. BONZANI, R. MONACO and M. G. ZAVATTARO Dipartimento di Matematica, Politecnico di Torino, 10129 Torino, Italia (Received August f987)
Communicated by E. Y. Rodin Abstract-This paper deals with the analysis of stochastic mathematical models in continuum mechanics. The first part of the paper provides a discussion on a general methodology for stochastic modelling of real systems in continuum physics and mechanics. The second part of the paper deals with models described by deterministic evolution partial differential equations with random initial boundary conditions in one space dimension and provides a mathematical method to compute the time evolution of the probability density in conjunction with the dependent variable.
1. A STOCHASTIC
MODEL
Consider a mathematical model dimension and in one dependent equation of the type
FOR
CLASSICAL
CONTINUUM
THEORY
of a system in continuum physics or mechanics in one space variable, such that its time-space evolution is described by an
au a34 $=f(t,x,u,~,~,~;r 1, at4
where r is a suitable new variable
set of parameters
into equation
characterizing
(l), one obtains
the model.
(1) Setting
a system of two first-order
au = u and casting
at
coupled
this
equations:
au (2)
The problem of defining some guidelines to obtain random models when the deterministic is given by equations (1) or (2) is dealt with in this section. The main sources of randomness in a continuum system can be listed as: (a) Random (b) Random (c) Random
limit
initial conditions: u,, = u(o, x; t = 0), u. = u(w, x; t = 0). boundary values: ubO = u(o, t; x = 0), ubl = u(w, t; x = 1). inner characteristics of the system and/or additional noise w(t, co).
Following the above scheme the stochastic system can be regarded as an input-output system, see Scheme 1 below, where initial and boundary conditions are the input, u, = u(w, x, t) is the output and the noise w,, if it exists, is also an input.
Scheme 1 207
I. BONZANI et al.
208
The randomness of the quantities defined in points (a) and (b) may be caused by the difficulty of observing and measuring them deterministically. A physically consistent origin of the nondeterministic measurability of the said quantities is indicated by the fact that if equations (2) are a mathematical model of a real physical system, then initial and boundary conditions will fluctuate (unless the observation time of the system itself is very short) about some mean values. One of the relevant reasons for these fluctuations lies in the fact that no real system can be regarded as fully isolated, so the outer world causes disturbances which need to be considered in terms of random variables [l]. In addition, as discussed in Refs [2] or [3], the medium of the continuum system can often be characterized by random properties. Bearing this in mind, we can consider a reasonable model for the initial and boundary conditions to be that defined by the Karhunen-Loeve expansion, say
(3) ub =
ubdt)
+ ~~v(d~&) Y
and analogously for the other quantities uo, ubo, ubl. In equations (3), M,(O) and by(o) are uncorrelated random variables with zero mean value, 4, and $, are known bounded functions of space and time, respectively, and aom and &, are the mean values of u. and ub, respectively. The number of terms v depends on the measures and estimates which we can realize for the quantities we are working on. Consider now the problem of the stochastic modelling of the system itself; i.e. of the term f in equations
(2). The randomness
can be localized
in the parameters
r and in the operators
au ah
-,
-.
ax ax2’
this randomness has to be regarded as a function of space, related to the structure of the matter characterizing the medium when the evolution defined by equation (1) occurs. In other words, we propose a mathematical model such that the characteristics of the medium are random (the physical justification for stochastic models of this kind is widely discussed in Refs [4-61). Consequently, each parameter ri of r can be modelled along the same lines as equations (3) as ri =
and the operator
au
-
ax
is replaced
+
CYidw)Pidx)
by the stochastic
operator
[
rim(X)
1 + Cu4PdX) Y
1
$
where yiy and 6, are random variables with zero mean value, py and py are known bounded functions of x and rim is the mean value of ri. Note that equation (5) is equivalent to assuming that the medium has an effective random length x*(o), such that
Then, the modelling
consists
of replacing
-=au
ax*
Similar
treatment
can be applied
g
by
au ax
ax ax*
N
a”[1
ax
+ a(o,x)].
to the second derivative,
when it appears
(7) in f.
A stochastic
model in continuum
mechanics
209
Finally, the noise w,, which simulates random disturbances of the outer world upon the system itself, can be related to a Wiener process, as is commonly done in the case of ordinary differential equations. This topic, however, is not developed further here. In fact, we only consider evolution equations with random parameters. Once the model has been defined, the output is the process u(w, x, t) solutions of the equations governing the system. Mathematical methods work towards a detailed calculation of u,. An important result to be studied consists of determining, for every x and t, the first and second probability densities P,(u, t, x)
and
P,(u”), u(‘), t,, t2; x).
In fact, if P, and P2 are known, all relevant statistical measures can be calculated. We now outline the contents of the forthcoming sections of this paper. Section 2 provides a mathematical description of a “relatively simpler” system belonging to the class defined in this section: a deterministic system with random initial boundary conditions. Section 3 presents the analysis of the problem; an approximated solution is obtained using an extension to the stochastic case [3, Chap. VI] of the Bellman differential quadrature method [7,Chap.X11]. Moreover, the approximated expressions for the first and the second probability densities related to the solution process are derived in Section 3. An example and a final discussion follow in Section 4.
2. MATHEMATICAL
FORMULATION
OF THE PROBLEM
Consider now, within the class of mathematical models described in the preceding section, the class of stochastic systems governed by a scalar differential equation with deterministic coefficients and random initial boundary conditions, defined on D = [0, I]. This class of problems, such that the randomness is limited to the initial boundary conditions, can be regarded as a relatively simpler class of problems with respect to the more general case where the randomness is localized in the parameters and operators also. Formally, it can be described by the equation
$=f(
au a5
t,x,u,~‘g
)
(8)
with constant boundary conditions and initial conditions: u(o,x,
t = 0) = u,(o,x);
(9)
u0 is a known bounded random variable defined in a complete probability space (Q g, P); R is the space of the elementary events o, g is a a-algebra of the subsets of $2 and P is the probability density. Remark
1
The random variable u. is assumed to be known in the m-order probability density P(u,(x,), for some fixed values xi, x2, . . . , x, of the space variable, where xi = 0 and . . .7 u,(x,)) x,= 1. The analysis will be carried out under the following hypothesis: uo(x2),
There are sufficient regularity conditions on f and on the initial boundary values, such that V(t,x) E [1 x D] and Vo ER the solution of equations (8) and (9) exists unique in the space
I. BONZANI et al.
210
@ is the space of the real functions of L,(R,W,P), respect to P, both with their first time derivative derivatives, with m 2 2.
continuous and bounded with and with their m-order space
The aim of this piece of research is to provide a mathematical method to derive, in approximated form, the time evolution of the first and second probability densities (see Ref. [S]) on the state variable u. The calculation of the probability densities provides a probabilistic time-space evolution of u, as discussed in Section 1. We point out that the analysis, when extended, can be applied to a wider class of equations than that introduced in equation (8), i.e. if (i) the dependent variable is a vector (u E [W”)and (ii) the system is modelled by an equation containing higher-order time derivatives. These extensions are rather technical and can be dealt with by generalizing the analysis developed in this paper-as discussed later.
3. ANALYSIS
OF THE
PROBLEM AND THE THE PROBABILITY
APPROXIMATED DENSITIES
DERIVATION
OF
Approximated expressions for the first and second probability densities related to the solution process of the physical system, mathematically described by equations (8) and (9) under the hypotheses defined in Section 2, will be derived here. The first step in reaching this goal is to determine the approximated solution process by means of the “stochastic interpolation method”, proposed in Refs [2; 3, Chap. VI] as a generalization to the stochastic case of the “Bellman differential quadrature method” [7, Chap. XII]. The application of the method consists of the following steps: (i) The space domain is discretized into a finite set of N points xi, with N < m. (ii) The solution u = u(o, t, x) satisfying equations (8) and (9) is written, assuming to be, for every PER, sufficiently smooth with respect to the space derivatives, as the interpolation
u(w~
t, x,
The pi(x), in equation Consequently,
N
t i=l
Pi(x)ui(07
Ui(W,
t);
(10) are suitable
the space derivatives
t)
=
u(w,
interpolation
appearing
in equation
t; x = Xi).
(10)
polynomials. (8) can be approximated
as follows:
and
(1lb) The coefficients aij, bij depend only on the choice of the polynomials pi(x), on the nodes N and on the measure of the spatial interval. If we choose, as in Ref. [2], the functions pi(x) as the (N - 1)th degree polynomials of Lagrangian type in the following form:
Pi(x)
where B’(xi) = dB(x)/dx
=
B(x) (x
and B(x) is a polynomial
-
Xi)B’(Xi)’
(12)
of degree N,
B(x) = (x - x1)(x - x2). . .(x - x,),
(13)
211
A stochastic model in continuum mechanics
the above coefficients
have the following
expressions:
i #j:
(14a)
and B”(Xi)
i =j:
ai.i= 2B’(x,)T
(14’4
where B”(x) = d2B(x)/dx2. Hence, under the assumption that equations (1 la, b) are valid, we have succeeded in reducing the original equation (8) into a system of (N - 2) ordinary differential equations on the random variable ui: i=2
,...,(N
- 1): 2 u =
with prescribed
initial
{“i>,
= fi(U, r; a, b), a
=
(15a) 4
{aij},
=
{b(j);
values Ui(O; t = 0) = UiO(W).
Equations (15a, b) in conjunction corresponding to i = 1 and i = N,
with the prescribed
u,(o;x
boundary
= 0) = u&D)
u,(o;x = 1) =
WW conditions
for x = 0 and x = 1,
\
(16)
u,,(w,,l
provide a system of N equations describing the time evolution in the N nodes of the random variable u. Obviously, the initial values in nodes 1 and N coincide with the boundary conditions for t = 0. More general boundary conditions on the space derivatives for u at the walls can be easily dealt with and result as prescribed linear combinations of Ui. This again reduces the number of independent ordinary differential equations to (N - 2). Remark
2
The analysis of problems in several dependent variable dimensions, outlined in point (i) of Section 2, can be carried out along the lines described in this section. In fact, in this case, u is a vector of R”, equation (8) is a system of n partial differential equations with n initial conditions and 2n boundary conditions. By applying the “stochastic interpolation method” the problem is reduced to consider n(N - 2) ordinary differential equations with prescribed initial boundary values. Remark
3
The problem corresponding to an equation containing higher-order time derivatives reduced to the preceeding one, as shown in equations (1) and (2) of Section 1. Remark
can be
4
Problems in the whole space XE R can be dealt with by resealing the space variable with a suitable non-linear transformation from the infinite space into the unit length strip. One needs in this case knowledge of the solution at infinity. Let us now deal with the evolution
equations
related to the random
initial conditions
(15a, b, 16).
212
1. BONZANI et al.
The analysis will be realized problem (15a, b, 16):
under
the following
hypotheses
on the solution
of the initial-value
The solution u(w, r) = @‘tu,(w) exists unique and continuous in uO, with continuous partial derivatives, a one-to-one mapping @‘tfor every realization of ~~(0). Accordingly, we can present random variables [S] :
the probability
density
using
Pl(U, t) = IJ(u -+ u,)lP,(u,
where J is the Jacobian
of the inverse
transformation
J(u + UJ =
The Jacobian J may be also computed equation [5,p. 1481, which leads to
IJI = exp
the classical
= a)-
and defines
formula
of change
‘u),
of
(17)
Q,; ‘, defined by
I2 I=
J(u,t).
(18)
by using a result derived
from the Liouville
evolution
1,
h being the r.h.s. of the set (15a), see also Ref. [3, Chap. IV] and Refs [8,9]. Equation (19) is particularly useful when the dimension of the set of the evolution equations is large; this is actually the case here because of the arbitrary choice of N. Furthermore, following Ref. [3, Chap. IV], we can write the second probability density P,, i.e. the probability density related to the random variables u(i) and u@), evaluated at times t = t, and t = t,, respectively:
P*(uCl),d2); t,, t2) = P,(u’l’; t1)[ 1 + 12(u”‘, II@);t,, t2)] - 1/Z
(20)
with
u(l) = u(t,; u&o)) = Qkl(o) 1,=11,
W=4
d2) = u(t,; u&o)) = @,u,(o) 11=,*
@lb)
I(u”‘, d2’; t,, t2) = J(lP’; t1)/J(d2’; t2).
(214
and
A stochastic model in continuum mechanics
213
Equation (21~) is the ratio between the Jacobians of the inverse map u -+ u0 calculated at t = t, and t = t,, respectively. Let us point out that equations (17) and (19) provide the first- and secondorder probability densities related to the N-dimensional variable u whose components are the approximated solutions of equation (8), computed in the N nodes. As a further step, one can calculate the probability densities related to each solution component ui. In fact, from the definition of marginal probability densities, one obtains, respectively:
Pl(Ui, t) =
Pl(U, t)du*
(22a)
s
and
P&&i, ~1~‘;t,, t2) =
W’b)
P,(u”‘, uc2); t,, t2)du’1)*du(2)*, s
where
du* = du,du2...dui_ldui+,...du,.
Of course one can use the solution of problem (15a, b, 16) in conjunction with interpolation (10) to compute, as in Ref. [2], the moments of the solution process. This is not the case for the analysis in this paper, which deals with the time evolution of the probability density. A practical application of this method is demonstrated in the problem considered in the next section.
4.
APPLICATION
TO
A MODEL
OF
THE
KINETIC
THEORY
OF
GASES
In this section a very simple application of the method developed above is proposed. application will be carried out essentially on tutorial grounds in order to show the practical of the procedure proposed. Consider the celebrated Broadwell [lo] model:
at4 au z
+
z
au
g
aw aw
_---_=v
at
ax
=
v2-
= uw -
This steps
uw,
v2,
2
-
uw,
(23)
which represents a dilute gas such that the particles have velocities with the same modulus and directions along the axes of an orthogonal frame 0(x, y, z). In particular, u(x, t) corresponds to the
I.
214
BONZANI et al.
numerical density of the particles with velocity (1, 0,0),w(x, t) to that with velocity (- l,O, 0), whereas u(x, t) is the numerical density for particles with velocities (0, l,O), (0, - l,O), (O,O, 1) or (O,O, - 1). Global existence and uniqueness for solutions of the initial-value problem was proved in Ref. [ 111. Let us attach to the set of equations (23) the random initial conditions
u(x, t = 0) = h + cr(w)sin nx = uO(x, CD), u(x, t = 0) = k, w(x, t = 0) = h,
(24)
where M(O) is a constant random variable related to a known PI(~). The problem (23,24) represents a gas in thermodynamical equilibrium subject to a random perturbation r(w)sin 7cx along the positive direction of the x-axis. By integrating equations (23) along the trajectories of the particles, starting from the initial conditions (24), after standard calculations, the problem (23,24) is reduced to the following single equation:
au au
sin[n(x - t)] + 2(h + Ic) U sin 71x I
- CkdX,~)-~l
g+x=
i
+
i
hl
Ca,(x,w) -
SinCdx- 4 + h + K* sin 7rx
with the initial condition given by the first of equations (24). Let us now apply to equation (25) the stochastic interpolation
gi_ &
-
1N
- j=
._{
C”idO)
, avuJ
-
hl
SiIl[7C(Xi
(25)
13
-
method.
We obtain,
t)]
sin71x, I
sinC4xi- Ql + h + K = Fi(u,t,oh sin rcxi
with initial
i = 1,. . . , N:
(26)
conditions
~(0)
= h + a(o)sinrrx,.
According to Remark 1, we assume that the probability density PI (uO) = P, (u,,, . . , uNo) is known and related to P, (a), i.e. the marginal probability P,(uio) must verify the condition
J
= PI (x)/sinzxi,
We can now calculate the probability equations (17, 19) we have
density
P,(u,t)
at every point
P,(u, t) = Mu + %)I P,(u,),
in time. In fact, applying
(27)
A stochastic
model in continuum
215
mechanics
where
1. After standard
calculations,
lJ(U
+
Uo)l
=
(28)
we obtain
exp
aiif + (UiO- h)
5 i
i=l
[
11
COS[7t(Xi- t)] + 2(h + IC)t sin nxi
.
(29)
of the Knowledge of the solution u(t; uo) and of the map u. = 0-i u allows explicit computation probability density P,(u, t). Further applications of equations (20)-(22) provide, finally, the expressions of the second probability density and of the marginal densities-as outlined at the end of the previous section. Before concluding this paper, let us consider a subexample of the one treated above which shows the rather good capability of the proposed method in applications to partial differential equations. In fact, consider the same model, equation (23), with the initial conditions
L&J =
h + a,
v. = k, w(J
=
h,
where tl is now a constant deterministic perturbation gas. In this case, equation (25) is reduced to
au au
of the thermodynamical
uo = h +
at+ax=A”+B,
equilibrium
of the
c7,
with
A=
-2(h+k+cr/2)
and
B=(h+k+01)2.
Considering the 1.h.s. of equation (30) as the total derivative du/dt of the numerical density of a gas whose particles move along the positive direction of the x-axis with unit velocity, then equation (30) itself has a unique positive exact solution (independent from x),
u = (h +
M +
B/A)exp(At)
- B/A,
for h2 - k2 + hct> 0. In Table 1 we show the time-dependent results for the exact values of u in comparison with those obtained by applying the stochastic interpolation method to equation (30) with N = 5, and integrating the corresponding set of ordinary differential equations with a standard fourth-order Runge-Kutta routine: the maximum error obtained is of the order of 10m8. Acknowledgements-This work has been carried out within the activities of the National Council for Research and supported by the Italian Ministry of Education and by the Strategical Project AITM of the GNFM.
(GNFM)
216
I. BONZANI et al.
Table I r=o x = 0 x = 0.25 x = 0.5 u = 0.75 .X=1
approx. approx. approx. approx. approx.
u u u u u
= = = = =
2 2 2 2 2
approx. approx. approx. ap!wx. approx.
u u u u u
= = = = =
1.65059711413 1.6505971141 1.65059711411 1.650597ll412 1.650597l1411
u u u u
= = = =
1.54535898157 1.54535898157 1.54535898155 1.54535898157 1.54535898159
Exact u = 2 I= x=0 Y= I = x = x=l
0.2 0.25 0.5 0.75
Exact u = 1.65059710596 f= x=0 x = x = Y= x =
0.4 0.25 0.5 0.75
approx. approx. approx. approx.
1
approx.u =
Exact u = 1.54535897665 t = 0.6 x=0 x = 0.25 x = 0.5 x = 0.75 x = 1
approx. approx. qqxox. approx. approx.
u u u u u
= = = = =
1.51366186341 1.51366186346 1.51366186345 1.51366186344 1.51366186345
u u u u
= = = =
1.50411487432 1.50411487441 1.50411487442 1.50411487442 1.50411487443
Exact u = 1.51366186122 f= x=0 x = x = x = x =
0.8 0.25 0.5 0.75
approx. approx. approx. approx.
1
approx.u =
Exact u = 1.50411487352
t=1 x=0
approx.u =
x = 0.25 x =0.5 x = 0.75 x=1
approx. approx. approx. approx.
u u u u
= = = =
1.50123937634 1.50123937639 1.50123937641 1.50123937641 1.50123937642
Exact u = 1.50123937609
REFERENCES 1. E. Nelson, Quanrum Fluctuations. Princeton Univ. Press, Princeton, N.J. (1985). 2. N. Bellomo, L. M. de Socio and R. Monaco, On the random heat equation: solutions by the stochastic adaptive method. Comput. Math. Applic. (1988). 3. N. Bellomo and R. Riganti, Nonlinear Stochastic Systems in Physics and Mechanics. World Scientific, Singapore (1987). 4. E. McShane, Stochastic Calculus and Stochastic Models. Academic Press, New York (1974). 5. T. T. Soong, Random Differential Equations in Science and Engineering. Academic Press, New York (1973). 6. G. Adomian, Stochastic Systems. Academic Press, New York (1983). 7. R. Bellman and G. Adomian, Partial Differential Equations. Reidel, Dordrecht (1985). 8. N. Bellomo and R. Riganti, Time evolution of the probability density and of the entropy function for a class of nonlinear systems in mathematical physics. Comput. Math. Applic. 12A, 663-676 (1986). 9. I. Bonzani, M. G. Zavattaro and N. Bellomo, On the continuous approximation of the probability density and of the entropy functions for nonlinear stochastic dynamical systems. Maths Comput. Simuln 29, (1987). 10. R. Gatignol, ThPorie Cinhtique des Gaz ci Repartition DiscrPte de Vitesses. Springer Lecture Notes in Physics, No. 36. Springer, Berlin (1975). 11. J. T. Beale, Large time behaviour of the Broadwell model of a discrete velocity gas. Commun. Math. Phys. 102, 217236 (1986).