- Email: [email protected]
A numerical scheme for the variance of the solution of the random transport equation
Applied Mathematics and Computation 190 (2007) 362–369 www.elsevier.com/locate/amc
A numerical scheme for the variance of the solution of the random transport equation Maria Cristina C. Cunha, Fa´bio Antonio Dorini
*
Departamento de Matema´tica Aplicada, IMECC – UNICAMP, Cidade Universita´ria Zeferino Vaz, CP. 6065, 13083-970 Campinas SP, Brazil
Abstract We present a numerical scheme, based on Godunov’s method (REA algorithm), for the variance of the solution of the 1D random linear transport equation, with homogeneous random velocity and stochastic initial condition. We obtain the stability conditions of the method and we also show its consistency with a deterministic nonhomogeneous advective–diffusive equation, which means convergency. Numerical results are considered to validate our scheme. Ó 2007 Elsevier Inc. All rights reserved. Keywords: Random linear transport equation; Finite volume schemes; Godunov’s method
1. Introduction In this work we are concerned about the variance of the solution of the random transport equation, Qt ðx; tÞ þ AQx ðx; tÞ ¼ 0; t > 0; x 2 R; Qðx; 0Þ ¼ Q0 ðxÞ;
ð1Þ
with a homogeneous random transport velocity, A, and stochastic initial condition, Q0(x). The solution, Q(x, t), is a random function. For the particular case, Riemann problem (1) with Q0 if x < 0; Qðx; 0Þ ¼ ð2Þ Qþ if x > 0; 0 þ where Q 0 and Q0 are random variables, we presented in [1] the expression for the solution: þ ð3Þ QR ðx; tÞ ¼ Q 0 þ X Q0 Q0 ; x x where X is a Bernoulli random variable with P ðX ¼ 0Þ ¼ 1 F A t and P ðX ¼ 1Þ ¼ F A t ; here FA(x) is the cumulative probability function of the random variable A.
*
Corresponding author. E-mail addresses: [email protected] (M.C.C. Cunha), [email protected] (F.A. Dorini).
0096-3003/$ - see front matter Ó 2007 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2007.01.038
M.C.C. Cunha, F.A. Dorini / Applied Mathematics and Computation 190 (2007) 362–369
363
þ Also, according to [1], considering the independence between A and both Q 0 , Q0 , the statistical mean and variance are given by x ð4Þ hQþ hQR ðx; tÞi ¼ hQ 0 i þ FA 0 i hQ0 i t and x xh xi 2 ð5Þ Var½QR ðx; tÞ ¼ Var½Q ½Var½Qþ 1 FA hQþ 0 þ FA 0 Var½Q0 þ F A 0 i hQ0 i : t t t In our point of view, the special attraction of (3)–(5) is their utilization in discretizations of stochastic equations, like (1). In [2] we present an explicit method to calculate the first statistical moment of Q(x, t), the solution of (1) with Qðx; 0Þ ¼ Q0 ðxÞ a random function. In that report we show that the Godunov method provides a numerical scheme for the statistical mean which is, under certain assumptions on the discretization, stable and consistent with a diffusive equation. Therefore, besides the scheme itself, the numerical approach also gives an effective equation compatible with one published in the literature. The aim of this paper is to improve the knowledge of the random solution of (1) with the random function Qðx; 0Þ ¼ Q0 ðxÞ. We present a numerical method to calculate the variance of Q(x, t), which is the quantity most commonly used to specify the dispersion of the distribution around its mean. In Section 2 we deduce the explicit numerical scheme using the Godunov’s ideas. Consistency, stability and convergency are analyzed in Section 3. Finally, in Section 4, we present some numerical examples.
2. The numerical scheme In this section we present the numerical scheme for the variance of the solution of (1). We denote the spatial and the time grid points by xj ¼ jDx and tn ¼ nDt, respectively, and the jth grid cell is Cj ¼ ðxj1=2 ; xjþ1=2 Þ, xj1=2 ¼ xj Dx . Let Qnj be an approximation of the cell average of Qðx; tn Þ: 2 Z Z xjþ1=2 1 1 n Qj ’ Qðx; tn Þdx ¼ Qðx; tn Þdx: ð6Þ Dx Cj Dx xj1=2 Assuming that the cell averages at time tn, Qnj , are known, we summarize the REA, for Reconstruct-EvolveAverage, algorithm [3,4] in three steps: e tn Þ, from the cell averages Qn . In our case we use [Step 1.] Reconstruct a piecewise polynomial function, Qðx; j n e tn Þ ¼ Qn for all x 2 Cj . the piecewise constant function with Qj in the jth cell, i.e., Qðx; j e tnþ1 Þ a time Dt [Step 2.] Evolve the equation exactly, or approximately, with this initial data to obtain Qðx; later. e tnþ1 Þ over each grid cell to obtain the new cell averages, i.e., [Step 3.] Average Qðx; Z 1 e tnþ1 Þdx: ¼ Qðx; Qnþ1 j Dx Cj At a time tn, the piecewise constant function, step 1, defines a set of Riemann problems in each x ¼ xj1=2 : the differential equation (1) with the initial condition ( n Qj1 if x < xj1=2 ; Qðx; tn Þ ¼ ð7Þ Qnj if x > xj1=2 : We may use (3) to find a local solution to each Riemann problem at a time
i x xj1=2 h n n Qðx; tnþ1=2 Þ ¼ Qj1 þ X Qj Qnj1 ; Dt=2
Dt 2
where, for a x sufficiently close to xj1=2 , X ðxÞ is the Bernoulli random variable: 1; P ðX ðxÞ ¼ 1Þ ¼ F A ðxÞ; X ðxÞ ¼ 0; P ðX ðxÞ ¼ 0Þ ¼ 1 F A ðxÞ:
later: ð8Þ
ð9Þ
364
M.C.C. Cunha, F.A. Dorini / Applied Mathematics and Computation 190 (2007) 362–369
Also, according with (4) and (5) and denoting Hj1=2 ðxÞ ¼ F A h i hQðx; tnþ1=2 Þi ¼ hQnj1 i þ Hj1=2 ðxÞ hQnj i hQnj1 i and
xxj1=2 Dt=2
, we have:
h i Var½Qðx; tnþ1=2 Þ ¼ Var½Qnj1 þ Hj1=2 ðxÞ Var½Qnj Var½Qnj1 h i2 þ Hj1=2 ðxÞð1 Hj1=2 ðxÞÞ hQnj i hQnj1 i :
ð10Þ
ð11Þ
e tnþ1=2 Þ, can be constructed by piecing together Therefore, the variance of the solution at time tnþ1=2 , Var½ Qðx; the local values of the variance, (11), provided that the half time step Dt2 is short enough such that adjacent Riemann problems do not interact between themselves. This requires that Dx and Dt must be chosen satisfying: Var½Qðxj1 ; tnþ1=2 Þ Var½Qnj1
and
Var½Qðxj ; tnþ1=2 Þ Var½Qnj ;
where the symbol ‘‘’’ means ‘‘sufficiently near to’’. Substituting these conditions in (11), the following conditions must be satisfied:
Dx FA 0 and Dt
FA
Dx 1: Dt
ð12Þ
; Dx Remark 1. We may regard (12) as a kind of CFL condition for the method. The interval Dx must Dt Dt Dx Dx contain the ‘‘effective support’’ of the density probability function of A. This means that outside Dt ; Dt the probability of A is sufficiently near to zero, i.e., it can be disregarded. The existence of an effective support is ensured by Chebyshev’s inequality: P fjA hAij P krA g 6 k12 , for all k > 0, where rA is the standard variation of A. Therefore, if we take k12 sufficiently close to zero, to escape from the interaction between solutions of Dt Riemann problems we must take ðjhAij þ krA Þ Dx 6 1. e tnþ1=2 Þ, x 2 ½xj1 ; xj ; its cell average will be Under the hypothesis (12), the expression (11) defines Var½ Qðx; denoted by Z xj 1 nþ1=2 e tnþ1=2 Þdx: V j1=2 ¼ Var½ Qðx; Dx xj1 Therefore, using (11) we have the cell average of the variance in ½xj1 ; xj at time t ¼ tnþ1=2 : Z xj n h io 1 nþ1=2 V nj1 þ Hj1=2 ðxÞ V nj V nj1 dx V j1=2 ¼ Dx xj1 Z xj h i2 1 Hj1=2 ðxÞð1 Hj1=2 ðxÞÞdx hQnj i hQnj1 i : þ Dx xj1 |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} C
Preliminary computational tests have shown that C reduces excessively the contribution of ½hQnj i hQnj1 i2 nþ1=2 to V j1=2 . The following approximation provides better results: Z xj 1 C¼ Hj1=2 ðxÞ½1 Hj1=2 ðxÞdx ’ Hj1=2 ðfÞ½1 Hj1=2 ðfÞ; Dx xj1 where f 2 ½xj1 ; xj is such that Hj1=2 ðfÞ½1 Hj1=2 ðfÞ ¼ max Hj1=2 ðxÞ½1 Hj1=2 ðxÞ: x2½xj1 ;xj
It is straightforward to show that f satisfies Hj1=2 ðfÞ ¼ 12. Thus, Hj1=2 ðfÞ½1 Hj1=2 ðfÞ ¼ 14 and, changing variables in the other integral, we have (Z Dx ) i2 h i 1h Dt Dt nþ1=2 n F A ðxÞdx V nj V nj1 þ hQnj i hQnj1 i : V j1=2 ¼ V j1 þ 2Dx 4 Dx Dt
ð13Þ
M.C.C. Cunha, F.A. Dorini / Applied Mathematics and Computation 190 (2007) 362–369
365
Lemma 2. Let A be a random variable and ½n; n an effective support of the density probability function, fA, of A, i.e., F A ðnÞ 0 and F A ðnÞ 1. Then Z n F A ðxÞdx n hAi: ð14Þ n
Proof. See [2]. h Dt Substituting (14) in (13) and denoting k ¼ Dx hAi, we have: i2 i kh i 1h 1h nþ1=2 V j1=2 ¼ V nj þ V nj1 V nj V nj1 þ hQnj i hQnj1 i : 2 2 4
ð15Þ
Now we can repeat the same procedure to obtain approximations of the solution in ½xj1=2 ; xjþ1=2 at the time tnþ1 : i2 i kh i 1h 1 h nþ1=2 nþ1=2 nþ1=2 nþ1=2 nþ1=2 n1=2 V jnþ1 ¼ V jþ1=2 þ V j1=2 V jþ1=2 V j1=2 þ hQjþ1=2 i hQj1=2 i : ð16Þ 2 2 4 The ideas of the Godunov method were also used in [2] to design a scheme for approximations of the statistical means of (1): i kh i 1h nþ1=2 hQj1=2 i ¼ hQnj i þ hQnj1 i hQnj i hQnj1 i ð17Þ 2 2 and hQnþ1 j i ¼
i kh i 1 h nþ1=2 nþ1=2 nþ1=2 nþ1=2 hQjþ1=2 i þ hQj1=2 i hQjþ1=2 i hQj1=2 i ; 2 2
or, joining these expressions, i 1 h i kh n n hQjþ1 i hQnj1 i þ ð1 þ k2 Þ hQnjþ1 i 2hQnj i þ hQnj1 i : hQnþ1 j i ¼ hQj i 2 4
ð18Þ
ð19Þ
Using (15) and (17) in (16), we can summarize the two step scheme for the variance in the explicit form: h i2 i 1 h i 1 kh V jnþ1 ¼ V nj V njþ1 V nj1 þ ð1 þ k2 Þ V njþ1 2V nj þ V nj1 þ ð1 kÞ hQnjþ1 i hQnj i 2 4 8 h i2 nh i h io2 1 1 hQnjþ1 i hQnj1 i k hQnjþ1 i 2hQnj i þ hQnj1 i þ ð1 þ kÞ hQnj i hQnj1 i þ ; ð20Þ 8 16 Dt hAi. where k ¼ Dx
3. Numerical analysis of the scheme In this section we analyze some numerical aspects of the method (19) and (20), for the mean and the variance of the solution of (1). We obtain the stability conditions of the scheme and we also show its consistency with a deterministic nonhomogeneous advective–diffusive system. 2
2 Proposition 3. For Dx Dt ¼ m fixed, the numerical scheme (19) and (20) yields an OðDx Þ approximation for uðx; tÞ and vðx; tÞ, solutions of the deterministic system of partial differential equations (PDE): ut þ hAiux ¼ 4m uxx ; ð21Þ vt þ hAivx ¼ 4m vxx þ 2m u2x :
Proof. Let vðx; tÞ and uðx; tÞ be smooth functions such that vðxj ; tn Þ ¼ V nj and uðxj ; tn Þ ¼ hQnj i. We have shown 2 in [2] that taking into account that Dx ¼ m is fixed, the numerical scheme (19) gives an OðDx2 Þ approximation Dt for uðx; tÞ, solution of the differential equation ut þ hAiux ¼ 4m uxx . Also, using the Taylor series in (20), we obtain:
366
M.C.C. Cunha, F.A. Dorini / Applied Mathematics and Computation 190 (2007) 362–369
Dt 2 vt þ vtt þ OðDt Þ þ hAi½vx þ OðDx2 Þ 2
2 m m Dx 2 2 2 uxx þ OðDx Þ ¼ ð1 þ k Þ½vxx þ OðDx Þ þ ð1 kÞ ux þ 4 8 2
2 k 2 m Dx m uxx þ OðDx2 Þ þ ux þ OðDx2 Þ uxx þ OðDx2 Þ þ ð1 þ kÞ ux : 8 2 4 2 2
Dt Since Dx ¼ m is fixed, we have k ¼ Dx hAi ¼ Dx hAi ¼ OðDxÞ and Dt ¼ OðDx2 Þ. Thus, grouping the terms of the m Dt same order, we have: m m vt þ hAivx ¼ vxx þ u2x þ OðDx2 Þ: 4 2
Remark 4. Although the consistency result is for any m, computational tests have shown that a good choice for 2 m ¼ Dx , in (21), to calculate hQðx; T Þi and VarðQðx; T ÞÞ, is m ¼ 2Var½AT . Therefore, the diffusive term, 4m, is well Dt approximated by 12 Var½AT . Remark 5. The modified equations in (21) constitute a decoupled deterministic nonhomogeneous convective– diffusive system whose transport terms are the mean of the velocity and the diffusive terms are the same. The source term in the second equation involves the spatial derivative of the mean, given by the first equation. Remark 6. In [2] we have shown that the stability condition of (19) is (12). On the other hand, since the terms corresponding to the mean can be considered source terms, the method for the variance, (20), has the same stability conditions, i.e., (12). As a linear problem, we have convergence. 4. Numerical examples To assess our method for the variance of the linear advective equation with random data we present two numerical examples. In Example 7 we solve a Riemann problem in which case the exact values of hQðx; tÞi and Var½Qðx; tÞ are known. In Example 8 we apply the method in a problem with random velocity and a correlated random field as the initial condition. In both examples we use A normally and lognormally distributed. The value of Dx is presented in the caption of the figures. The value of Dt was chosen based on Remark 4, i.e., we 2 used m ¼ Dx ¼ 2Var½AT . All the numerical experiments presented in this section were computed in double preDt cision with some MATLAB codes on a 3.0 GHz Pentium 4 with 512 MB of memory. Example 7. Let us consider the random PDE (1) with the mean and the variance of the initial condition given by
Mean of solution
Variance of solution Proposed method Exact solution
1.2
0.5
1
0.45
0.8
0.4
0.6
0.35
0.4
0.3
0.2
0.25
0
0.2
–0.2
0.15 –1
0
1
2
3
–1
Fig. 1. Dx = 0.02 and T = 0.3.
0
1
2
3
M.C.C. Cunha, F.A. Dorini / Applied Mathematics and Computation 190 (2007) 362–369
hQðx; 0Þi ¼
1
if x < 0
0
if x P 0
Var½Qðx; 0Þ ¼
and
0:16
if x < 0;
0:25
if x P 0:
367
In Figs. 1 and 4 we compare the approximations of the mean and the variance calculated using (19) and (20), respectively, with the exact values given by (4) and (5). We plot the results in T = 0.3 and T = 0.5. To observe the influence of the velocity variation we use two models: (i) A normally distributed, A = N(1.0, 0.6), in Figs. 1 and 2; (ii) A lognormally distributed, A ¼ expðnÞ, n = N(0.5, 0.25), in Figs. 3 and 4.
Mean of solution
Variance of solution Proposed method Exact solution
1.2
0.5
1
0.45
0.8
0.4
0.6
0.35
0.4
0.3
0.2
0.25 0.2
0 –0.2
0.15 –1
0
1
2
3
–1
0
1
2
3
Fig. 2. Dx = 0.02 and T = 0.5.
Mean of solution
Variance of solution Proposed method Exact solution
1.2 1
0.5 0.45
0.8
0.4
0.6
0.35
0.4
0.3
0.2
0.25 0.2
0 –0.2
0.15 –1
0
1
2
3
–1
0
1
2
3
Fig. 3. Dx = 0.01 and T = 0.3.
Mean of solution
Variance of solution Proposed method Exact solution
1.2
0.5
1
0.45
0.8
0.4
0.6
0.35
0.4
0.3
0.2
0.25
0
0.2
–0.2
0.15 –1
0
1
2
3
–1
Fig. 4. Dx = 0.01 and T = 0.5.
0
1
2
3
368
M.C.C. Cunha, F.A. Dorini / Applied Mathematics and Computation 190 (2007) 362–369
Example 8. In this example we take the random PDE (1) with initial condition, Q0(x), being the random field with mean ( 1; x 2 ð1:4; 2:2Þ; hQ0 ðxÞi ¼ ð22Þ 20ðx0:25Þ2 e ; otherwise Mean of solution
Mean of initial condition Proposed method Monte Carlo simulations
1.2
Variance of solution
0.4
1 0.35 0.8 0.3
0.6
0.25
0.4
0.2
0.2
0.15
0 –0.2 –1
0.1 –1
3
2
1
0
3
2
1
0
Fig. 5. Dx = 0.02 and T = 0.3.
Mean of solution
Mean of initial condition Proposed method Monte Carlo simulations
1
Variance of solution 0.45 0.4
0.8 0.35 0.6
0.3
0.4
0.25
0.2
0.2
0
0.15
–0.2 –1
0
1
2
0.1 –1
3
0
1
2
3
Fig. 6. Dx = 0.02 and T = 0.5.
Mean of solution
Mean of initial condition Proposed method Monte Carlo simulations
1.2
Variance of solution 0.45 0.4
1 0.35
0.8
0.3
0.6
0.25
0.4
0.2
0.2
0.15
0 –0.2 –0.5
0.1 0
0.5
1
1.5
2
2.5
3
3.5
–0.5
Fig. 7. Dx = 0.01 and T = 0.3.
0
0.5
1
1.5
2
2.5
3
3.5
M.C.C. Cunha, F.A. Dorini / Applied Mathematics and Computation 190 (2007) 362–369 Mean of solution
Mean of initial condition Proposed method Monte Carlo simulations
1.2
Variance of solution 0.45 0.4
1
0.35
0.8
0.3
0.6
0.25
0.4
0.2
0.2
0.15
0 –0.2 –0.5
369
0.1 0
0.5
1
1.5
2
2.5
3
3.5
–0.5
0
0.5
1
1.5
2
2.5
3
3.5
Fig. 8. Dx = 0.01 and T = 0.5.
and covariance Covðx; ~xÞ ¼ r2 expðbjx ~xjÞ, with Var½Q0 ðxÞ ¼ r2 constant; the parameter b > 0 governs the decay rate of the spatial correlation. In our tests we use b = 40 and r2 = 0.12. The numerical results are compared with Monte Carlo simulations using suites of realizations of A and Q0(x), with A and Q0(x) independents. As it is known, the analytical solution of each realization A(x) and Q0(x, x) is given by Qðx; t; xÞ ¼ Q0 ðx AðxÞt; xÞ. The 2000 realizations of the correlated random field Q0(x) are generated using the matrix decomposition method, a direct method for generating correlated random fields (for example [5], Chapter 3). As in the previous example we use two models of velocity: (i) A normally distributed, A = N(0.5, 0.6), in Figs. 5 and 6; (ii) A lognormally distributed, A = exp(n), n = N(0.15, 0.25), in Figs. 7 and 8. 5. Concluding remarks In this paper we extend the ideas presented in our previous work [2] to obtain more information about the statistical properties of the solution to one dimensional stochastic transport partial differential equations. We show that the ideas of the Godunov method can also be used to design a numerical scheme to calculate the variance of the solution: (19) and (20). We also present the stability conditions and the consistency of the numerical scheme with the decoupled system of convective–diffusive equations (21). Computational results are confronted with the exact solution, in the Riemann problem, and with Monte Carlo simulations in a more general situation. As far as we know, this kind of methodology has not been used to deal with differential equations with uncertainties in the parameters. This approach can represent a gain if compared with the Monte Carlo simulations, the effective equations, or other usual methodologies. Extensions to more general situations, for example the variable velocity case and 2D problems will be presented in the future. Acknowledgement Our acknowledgments to the Brazilian Council for Development of Science and Technology (CNPq) through the Grant 5551463/02-3 and doctoral scholarship. References [1] M.C.C. Cunha, F.A. Dorini, Riemann problem for random advective equation. Technical Report 17/2006, State University of Campinas, Department of Applied Mathematics, Institute of Mathematics, Campinas, SP, 2006.. [2] F.A. Dorini, M.C.C. Cunha, A finite volume method for the mean of the solution of the random transport equation, Appl. Math. Comput. (2006), doi:10.1016/j.amc.2006.09.029. [3] S.K. Godunov, A difference method for numerical calculation of discontinuous solutions of the equations of hydrodynamics, Mat. Sb. 47 (1959) 271–306. [4] R.J. LeVeque, Finite Volume Methods for Hyperbolic Problems, Cambridge University Press, Cambridge, 2002. [5] D. Zhang, Stochastic Methods for Flows in Porous Media – Coping with Uncertainties, Academic Press, 2002.
A numerical scheme for the variance of the solution of the random transport equation Maria Cristina C. Cunha, Fa´bio Antonio Dorini
*
Departamento de Matema´tica Aplicada, IMECC – UNICAMP, Cidade Universita´ria Zeferino Vaz, CP. 6065, 13083-970 Campinas SP, Brazil
Abstract We present a numerical scheme, based on Godunov’s method (REA algorithm), for the variance of the solution of the 1D random linear transport equation, with homogeneous random velocity and stochastic initial condition. We obtain the stability conditions of the method and we also show its consistency with a deterministic nonhomogeneous advective–diffusive equation, which means convergency. Numerical results are considered to validate our scheme. Ó 2007 Elsevier Inc. All rights reserved. Keywords: Random linear transport equation; Finite volume schemes; Godunov’s method
1. Introduction In this work we are concerned about the variance of the solution of the random transport equation, Qt ðx; tÞ þ AQx ðx; tÞ ¼ 0; t > 0; x 2 R; Qðx; 0Þ ¼ Q0 ðxÞ;
ð1Þ
with a homogeneous random transport velocity, A, and stochastic initial condition, Q0(x). The solution, Q(x, t), is a random function. For the particular case, Riemann problem (1) with Q0 if x < 0; Qðx; 0Þ ¼ ð2Þ Qþ if x > 0; 0 þ where Q 0 and Q0 are random variables, we presented in [1] the expression for the solution: þ ð3Þ QR ðx; tÞ ¼ Q 0 þ X Q0 Q0 ; x x where X is a Bernoulli random variable with P ðX ¼ 0Þ ¼ 1 F A t and P ðX ¼ 1Þ ¼ F A t ; here FA(x) is the cumulative probability function of the random variable A.
*
Corresponding author. E-mail addresses: [email protected] (M.C.C. Cunha), [email protected] (F.A. Dorini).
0096-3003/$ - see front matter Ó 2007 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2007.01.038
M.C.C. Cunha, F.A. Dorini / Applied Mathematics and Computation 190 (2007) 362–369
363
þ Also, according to [1], considering the independence between A and both Q 0 , Q0 , the statistical mean and variance are given by x ð4Þ hQþ hQR ðx; tÞi ¼ hQ 0 i þ FA 0 i hQ0 i t and x xh xi 2 ð5Þ Var½QR ðx; tÞ ¼ Var½Q ½Var½Qþ 1 FA hQþ 0 þ FA 0 Var½Q0 þ F A 0 i hQ0 i : t t t In our point of view, the special attraction of (3)–(5) is their utilization in discretizations of stochastic equations, like (1). In [2] we present an explicit method to calculate the first statistical moment of Q(x, t), the solution of (1) with Qðx; 0Þ ¼ Q0 ðxÞ a random function. In that report we show that the Godunov method provides a numerical scheme for the statistical mean which is, under certain assumptions on the discretization, stable and consistent with a diffusive equation. Therefore, besides the scheme itself, the numerical approach also gives an effective equation compatible with one published in the literature. The aim of this paper is to improve the knowledge of the random solution of (1) with the random function Qðx; 0Þ ¼ Q0 ðxÞ. We present a numerical method to calculate the variance of Q(x, t), which is the quantity most commonly used to specify the dispersion of the distribution around its mean. In Section 2 we deduce the explicit numerical scheme using the Godunov’s ideas. Consistency, stability and convergency are analyzed in Section 3. Finally, in Section 4, we present some numerical examples.
2. The numerical scheme In this section we present the numerical scheme for the variance of the solution of (1). We denote the spatial and the time grid points by xj ¼ jDx and tn ¼ nDt, respectively, and the jth grid cell is Cj ¼ ðxj1=2 ; xjþ1=2 Þ, xj1=2 ¼ xj Dx . Let Qnj be an approximation of the cell average of Qðx; tn Þ: 2 Z Z xjþ1=2 1 1 n Qj ’ Qðx; tn Þdx ¼ Qðx; tn Þdx: ð6Þ Dx Cj Dx xj1=2 Assuming that the cell averages at time tn, Qnj , are known, we summarize the REA, for Reconstruct-EvolveAverage, algorithm [3,4] in three steps: e tn Þ, from the cell averages Qn . In our case we use [Step 1.] Reconstruct a piecewise polynomial function, Qðx; j n e tn Þ ¼ Qn for all x 2 Cj . the piecewise constant function with Qj in the jth cell, i.e., Qðx; j e tnþ1 Þ a time Dt [Step 2.] Evolve the equation exactly, or approximately, with this initial data to obtain Qðx; later. e tnþ1 Þ over each grid cell to obtain the new cell averages, i.e., [Step 3.] Average Qðx; Z 1 e tnþ1 Þdx: ¼ Qðx; Qnþ1 j Dx Cj At a time tn, the piecewise constant function, step 1, defines a set of Riemann problems in each x ¼ xj1=2 : the differential equation (1) with the initial condition ( n Qj1 if x < xj1=2 ; Qðx; tn Þ ¼ ð7Þ Qnj if x > xj1=2 : We may use (3) to find a local solution to each Riemann problem at a time
i x xj1=2 h n n Qðx; tnþ1=2 Þ ¼ Qj1 þ X Qj Qnj1 ; Dt=2
Dt 2
where, for a x sufficiently close to xj1=2 , X ðxÞ is the Bernoulli random variable: 1; P ðX ðxÞ ¼ 1Þ ¼ F A ðxÞ; X ðxÞ ¼ 0; P ðX ðxÞ ¼ 0Þ ¼ 1 F A ðxÞ:
later: ð8Þ
ð9Þ
364
M.C.C. Cunha, F.A. Dorini / Applied Mathematics and Computation 190 (2007) 362–369
Also, according with (4) and (5) and denoting Hj1=2 ðxÞ ¼ F A h i hQðx; tnþ1=2 Þi ¼ hQnj1 i þ Hj1=2 ðxÞ hQnj i hQnj1 i and
xxj1=2 Dt=2
, we have:
h i Var½Qðx; tnþ1=2 Þ ¼ Var½Qnj1 þ Hj1=2 ðxÞ Var½Qnj Var½Qnj1 h i2 þ Hj1=2 ðxÞð1 Hj1=2 ðxÞÞ hQnj i hQnj1 i :
ð10Þ
ð11Þ
e tnþ1=2 Þ, can be constructed by piecing together Therefore, the variance of the solution at time tnþ1=2 , Var½ Qðx; the local values of the variance, (11), provided that the half time step Dt2 is short enough such that adjacent Riemann problems do not interact between themselves. This requires that Dx and Dt must be chosen satisfying: Var½Qðxj1 ; tnþ1=2 Þ Var½Qnj1
and
Var½Qðxj ; tnþ1=2 Þ Var½Qnj ;
where the symbol ‘‘’’ means ‘‘sufficiently near to’’. Substituting these conditions in (11), the following conditions must be satisfied:
Dx FA 0 and Dt
FA
Dx 1: Dt
ð12Þ
; Dx Remark 1. We may regard (12) as a kind of CFL condition for the method. The interval Dx must Dt Dt Dx Dx contain the ‘‘effective support’’ of the density probability function of A. This means that outside Dt ; Dt the probability of A is sufficiently near to zero, i.e., it can be disregarded. The existence of an effective support is ensured by Chebyshev’s inequality: P fjA hAij P krA g 6 k12 , for all k > 0, where rA is the standard variation of A. Therefore, if we take k12 sufficiently close to zero, to escape from the interaction between solutions of Dt Riemann problems we must take ðjhAij þ krA Þ Dx 6 1. e tnþ1=2 Þ, x 2 ½xj1 ; xj ; its cell average will be Under the hypothesis (12), the expression (11) defines Var½ Qðx; denoted by Z xj 1 nþ1=2 e tnþ1=2 Þdx: V j1=2 ¼ Var½ Qðx; Dx xj1 Therefore, using (11) we have the cell average of the variance in ½xj1 ; xj at time t ¼ tnþ1=2 : Z xj n h io 1 nþ1=2 V nj1 þ Hj1=2 ðxÞ V nj V nj1 dx V j1=2 ¼ Dx xj1 Z xj h i2 1 Hj1=2 ðxÞð1 Hj1=2 ðxÞÞdx hQnj i hQnj1 i : þ Dx xj1 |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} C
Preliminary computational tests have shown that C reduces excessively the contribution of ½hQnj i hQnj1 i2 nþ1=2 to V j1=2 . The following approximation provides better results: Z xj 1 C¼ Hj1=2 ðxÞ½1 Hj1=2 ðxÞdx ’ Hj1=2 ðfÞ½1 Hj1=2 ðfÞ; Dx xj1 where f 2 ½xj1 ; xj is such that Hj1=2 ðfÞ½1 Hj1=2 ðfÞ ¼ max Hj1=2 ðxÞ½1 Hj1=2 ðxÞ: x2½xj1 ;xj
It is straightforward to show that f satisfies Hj1=2 ðfÞ ¼ 12. Thus, Hj1=2 ðfÞ½1 Hj1=2 ðfÞ ¼ 14 and, changing variables in the other integral, we have (Z Dx ) i2 h i 1h Dt Dt nþ1=2 n F A ðxÞdx V nj V nj1 þ hQnj i hQnj1 i : V j1=2 ¼ V j1 þ 2Dx 4 Dx Dt
ð13Þ
M.C.C. Cunha, F.A. Dorini / Applied Mathematics and Computation 190 (2007) 362–369
365
Lemma 2. Let A be a random variable and ½n; n an effective support of the density probability function, fA, of A, i.e., F A ðnÞ 0 and F A ðnÞ 1. Then Z n F A ðxÞdx n hAi: ð14Þ n
Proof. See [2]. h Dt Substituting (14) in (13) and denoting k ¼ Dx hAi, we have: i2 i kh i 1h 1h nþ1=2 V j1=2 ¼ V nj þ V nj1 V nj V nj1 þ hQnj i hQnj1 i : 2 2 4
ð15Þ
Now we can repeat the same procedure to obtain approximations of the solution in ½xj1=2 ; xjþ1=2 at the time tnþ1 : i2 i kh i 1h 1 h nþ1=2 nþ1=2 nþ1=2 nþ1=2 nþ1=2 n1=2 V jnþ1 ¼ V jþ1=2 þ V j1=2 V jþ1=2 V j1=2 þ hQjþ1=2 i hQj1=2 i : ð16Þ 2 2 4 The ideas of the Godunov method were also used in [2] to design a scheme for approximations of the statistical means of (1): i kh i 1h nþ1=2 hQj1=2 i ¼ hQnj i þ hQnj1 i hQnj i hQnj1 i ð17Þ 2 2 and hQnþ1 j i ¼
i kh i 1 h nþ1=2 nþ1=2 nþ1=2 nþ1=2 hQjþ1=2 i þ hQj1=2 i hQjþ1=2 i hQj1=2 i ; 2 2
or, joining these expressions, i 1 h i kh n n hQjþ1 i hQnj1 i þ ð1 þ k2 Þ hQnjþ1 i 2hQnj i þ hQnj1 i : hQnþ1 j i ¼ hQj i 2 4
ð18Þ
ð19Þ
Using (15) and (17) in (16), we can summarize the two step scheme for the variance in the explicit form: h i2 i 1 h i 1 kh V jnþ1 ¼ V nj V njþ1 V nj1 þ ð1 þ k2 Þ V njþ1 2V nj þ V nj1 þ ð1 kÞ hQnjþ1 i hQnj i 2 4 8 h i2 nh i h io2 1 1 hQnjþ1 i hQnj1 i k hQnjþ1 i 2hQnj i þ hQnj1 i þ ð1 þ kÞ hQnj i hQnj1 i þ ; ð20Þ 8 16 Dt hAi. where k ¼ Dx
3. Numerical analysis of the scheme In this section we analyze some numerical aspects of the method (19) and (20), for the mean and the variance of the solution of (1). We obtain the stability conditions of the scheme and we also show its consistency with a deterministic nonhomogeneous advective–diffusive system. 2
2 Proposition 3. For Dx Dt ¼ m fixed, the numerical scheme (19) and (20) yields an OðDx Þ approximation for uðx; tÞ and vðx; tÞ, solutions of the deterministic system of partial differential equations (PDE): ut þ hAiux ¼ 4m uxx ; ð21Þ vt þ hAivx ¼ 4m vxx þ 2m u2x :
Proof. Let vðx; tÞ and uðx; tÞ be smooth functions such that vðxj ; tn Þ ¼ V nj and uðxj ; tn Þ ¼ hQnj i. We have shown 2 in [2] that taking into account that Dx ¼ m is fixed, the numerical scheme (19) gives an OðDx2 Þ approximation Dt for uðx; tÞ, solution of the differential equation ut þ hAiux ¼ 4m uxx . Also, using the Taylor series in (20), we obtain:
366
M.C.C. Cunha, F.A. Dorini / Applied Mathematics and Computation 190 (2007) 362–369
Dt 2 vt þ vtt þ OðDt Þ þ hAi½vx þ OðDx2 Þ 2
2 m m Dx 2 2 2 uxx þ OðDx Þ ¼ ð1 þ k Þ½vxx þ OðDx Þ þ ð1 kÞ ux þ 4 8 2
2 k 2 m Dx m uxx þ OðDx2 Þ þ ux þ OðDx2 Þ uxx þ OðDx2 Þ þ ð1 þ kÞ ux : 8 2 4 2 2
Dt Since Dx ¼ m is fixed, we have k ¼ Dx hAi ¼ Dx hAi ¼ OðDxÞ and Dt ¼ OðDx2 Þ. Thus, grouping the terms of the m Dt same order, we have: m m vt þ hAivx ¼ vxx þ u2x þ OðDx2 Þ: 4 2
Remark 4. Although the consistency result is for any m, computational tests have shown that a good choice for 2 m ¼ Dx , in (21), to calculate hQðx; T Þi and VarðQðx; T ÞÞ, is m ¼ 2Var½AT . Therefore, the diffusive term, 4m, is well Dt approximated by 12 Var½AT . Remark 5. The modified equations in (21) constitute a decoupled deterministic nonhomogeneous convective– diffusive system whose transport terms are the mean of the velocity and the diffusive terms are the same. The source term in the second equation involves the spatial derivative of the mean, given by the first equation. Remark 6. In [2] we have shown that the stability condition of (19) is (12). On the other hand, since the terms corresponding to the mean can be considered source terms, the method for the variance, (20), has the same stability conditions, i.e., (12). As a linear problem, we have convergence. 4. Numerical examples To assess our method for the variance of the linear advective equation with random data we present two numerical examples. In Example 7 we solve a Riemann problem in which case the exact values of hQðx; tÞi and Var½Qðx; tÞ are known. In Example 8 we apply the method in a problem with random velocity and a correlated random field as the initial condition. In both examples we use A normally and lognormally distributed. The value of Dx is presented in the caption of the figures. The value of Dt was chosen based on Remark 4, i.e., we 2 used m ¼ Dx ¼ 2Var½AT . All the numerical experiments presented in this section were computed in double preDt cision with some MATLAB codes on a 3.0 GHz Pentium 4 with 512 MB of memory. Example 7. Let us consider the random PDE (1) with the mean and the variance of the initial condition given by
Mean of solution
Variance of solution Proposed method Exact solution
1.2
0.5
1
0.45
0.8
0.4
0.6
0.35
0.4
0.3
0.2
0.25
0
0.2
–0.2
0.15 –1
0
1
2
3
–1
Fig. 1. Dx = 0.02 and T = 0.3.
0
1
2
3
M.C.C. Cunha, F.A. Dorini / Applied Mathematics and Computation 190 (2007) 362–369
hQðx; 0Þi ¼
1
if x < 0
0
if x P 0
Var½Qðx; 0Þ ¼
and
0:16
if x < 0;
0:25
if x P 0:
367
In Figs. 1 and 4 we compare the approximations of the mean and the variance calculated using (19) and (20), respectively, with the exact values given by (4) and (5). We plot the results in T = 0.3 and T = 0.5. To observe the influence of the velocity variation we use two models: (i) A normally distributed, A = N(1.0, 0.6), in Figs. 1 and 2; (ii) A lognormally distributed, A ¼ expðnÞ, n = N(0.5, 0.25), in Figs. 3 and 4.
Mean of solution
Variance of solution Proposed method Exact solution
1.2
0.5
1
0.45
0.8
0.4
0.6
0.35
0.4
0.3
0.2
0.25 0.2
0 –0.2
0.15 –1
0
1
2
3
–1
0
1
2
3
Fig. 2. Dx = 0.02 and T = 0.5.
Mean of solution
Variance of solution Proposed method Exact solution
1.2 1
0.5 0.45
0.8
0.4
0.6
0.35
0.4
0.3
0.2
0.25 0.2
0 –0.2
0.15 –1
0
1
2
3
–1
0
1
2
3
Fig. 3. Dx = 0.01 and T = 0.3.
Mean of solution
Variance of solution Proposed method Exact solution
1.2
0.5
1
0.45
0.8
0.4
0.6
0.35
0.4
0.3
0.2
0.25
0
0.2
–0.2
0.15 –1
0
1
2
3
–1
Fig. 4. Dx = 0.01 and T = 0.5.
0
1
2
3
368
M.C.C. Cunha, F.A. Dorini / Applied Mathematics and Computation 190 (2007) 362–369
Example 8. In this example we take the random PDE (1) with initial condition, Q0(x), being the random field with mean ( 1; x 2 ð1:4; 2:2Þ; hQ0 ðxÞi ¼ ð22Þ 20ðx0:25Þ2 e ; otherwise Mean of solution
Mean of initial condition Proposed method Monte Carlo simulations
1.2
Variance of solution
0.4
1 0.35 0.8 0.3
0.6
0.25
0.4
0.2
0.2
0.15
0 –0.2 –1
0.1 –1
3
2
1
0
3
2
1
0
Fig. 5. Dx = 0.02 and T = 0.3.
Mean of solution
Mean of initial condition Proposed method Monte Carlo simulations
1
Variance of solution 0.45 0.4
0.8 0.35 0.6
0.3
0.4
0.25
0.2
0.2
0
0.15
–0.2 –1
0
1
2
0.1 –1
3
0
1
2
3
Fig. 6. Dx = 0.02 and T = 0.5.
Mean of solution
Mean of initial condition Proposed method Monte Carlo simulations
1.2
Variance of solution 0.45 0.4
1 0.35
0.8
0.3
0.6
0.25
0.4
0.2
0.2
0.15
0 –0.2 –0.5
0.1 0
0.5
1
1.5
2
2.5
3
3.5
–0.5
Fig. 7. Dx = 0.01 and T = 0.3.
0
0.5
1
1.5
2
2.5
3
3.5
M.C.C. Cunha, F.A. Dorini / Applied Mathematics and Computation 190 (2007) 362–369 Mean of solution
Mean of initial condition Proposed method Monte Carlo simulations
1.2
Variance of solution 0.45 0.4
1
0.35
0.8
0.3
0.6
0.25
0.4
0.2
0.2
0.15
0 –0.2 –0.5
369
0.1 0
0.5
1
1.5
2
2.5
3
3.5
–0.5
0
0.5
1
1.5
2
2.5
3
3.5
Fig. 8. Dx = 0.01 and T = 0.5.
and covariance Covðx; ~xÞ ¼ r2 expðbjx ~xjÞ, with Var½Q0 ðxÞ ¼ r2 constant; the parameter b > 0 governs the decay rate of the spatial correlation. In our tests we use b = 40 and r2 = 0.12. The numerical results are compared with Monte Carlo simulations using suites of realizations of A and Q0(x), with A and Q0(x) independents. As it is known, the analytical solution of each realization A(x) and Q0(x, x) is given by Qðx; t; xÞ ¼ Q0 ðx AðxÞt; xÞ. The 2000 realizations of the correlated random field Q0(x) are generated using the matrix decomposition method, a direct method for generating correlated random fields (for example [5], Chapter 3). As in the previous example we use two models of velocity: (i) A normally distributed, A = N(0.5, 0.6), in Figs. 5 and 6; (ii) A lognormally distributed, A = exp(n), n = N(0.15, 0.25), in Figs. 7 and 8. 5. Concluding remarks In this paper we extend the ideas presented in our previous work [2] to obtain more information about the statistical properties of the solution to one dimensional stochastic transport partial differential equations. We show that the ideas of the Godunov method can also be used to design a numerical scheme to calculate the variance of the solution: (19) and (20). We also present the stability conditions and the consistency of the numerical scheme with the decoupled system of convective–diffusive equations (21). Computational results are confronted with the exact solution, in the Riemann problem, and with Monte Carlo simulations in a more general situation. As far as we know, this kind of methodology has not been used to deal with differential equations with uncertainties in the parameters. This approach can represent a gain if compared with the Monte Carlo simulations, the effective equations, or other usual methodologies. Extensions to more general situations, for example the variable velocity case and 2D problems will be presented in the future. Acknowledgement Our acknowledgments to the Brazilian Council for Development of Science and Technology (CNPq) through the Grant 5551463/02-3 and doctoral scholarship. References [1] M.C.C. Cunha, F.A. Dorini, Riemann problem for random advective equation. Technical Report 17/2006, State University of Campinas, Department of Applied Mathematics, Institute of Mathematics, Campinas, SP, 2006.