Skip to main content
Log in

A multilevel Monte Carlo method for asymptotic-preserving particle schemes in the diffusive limit

  • Published:
Numerische Mathematik Aims and scope Submit manuscript

Abstract

Kinetic equations model distributions of particles in position-velocity phase space. Often, one is interested in studying the long-time behavior of particles in high-collisional regimes in which an approximate (advection)-diffusion model holds. In this paper we consider the diffusive scaling. Classical particle-based techniques suffer from a strict time-step restriction in this limit, to maintain stability. Asymptotic-preserving schemes avoid this problem, but introduce an additional time discretization error, possibly resulting in an unacceptably large bias for larger time steps. Here, we present and analyze a multilevel Monte Carlo scheme that reduces this bias by combining estimates using a hierarchy of different time step sizes. We demonstrate how to correlate trajectories from this scheme, using different time steps. We also present a strategy for selecting the levels in the multilevel scheme. Our approach significantly reduces the computation required to perform accurate simulations of the considered kinetic equations, compared to classical Monte Carlo approaches.

This is a preview of subscription content, log in via an institution to check access.

Access this article

Price excludes VAT (USA)
Tax calculation will be finalised during checkout.

Instant access to the full article PDF.

Institutional subscriptions

Fig. 1
Fig. 2
Fig. 3
Fig. 4
Fig. 5
Fig. 6
Fig. 7
Fig. 8
Fig. 9

Similar content being viewed by others

References

  1. Anderson, D.F., Higham, D.J.: Multilevel Monte Carlo for continuous time Markov chains, with applications in biochemical kinetics. Multiscale Model. Simul. 10(1), 146–179 (2012)

    Article  MathSciNet  Google Scholar 

  2. Bennoune, M., Lemou, M., Mieussens, L.: Uniformly stable numerical schemes for the Boltzmann equation preserving the compressible Navier–Stokes asymptotics. J. Comput. Phys. 227(8), 3781–3803 (2008)

    Article  MathSciNet  Google Scholar 

  3. Bhatnagar, P.L., Gross, E.P., Krook, M.: A model for collision processes in gases. I. Small amplitude processes in charged and neutral one-component systems. Phys. Rev. 94(3), 511–525 (1954)

    Article  Google Scholar 

  4. Birdsall, C.K., Langdon, A.B.: Plasma Physics via Computer Simulation. Series in Plasma Physics and Fluid Dynamics. Taylor & Francis, New York (2004)

    Google Scholar 

  5. Blondeel, P., Robbe, P., Van Hoorickx, C., François, S., Lombaert, G., Vandewalle, S.: p-refined multilevel quasi-Monte Carlo for Galerkin finite element methods with applications in civil engineering. Algorithms 13(5), 110 (2020)

    Article  MathSciNet  Google Scholar 

  6. Boscarino, S., Pareschi, L., Russo, G.: Implicit–explicit Runge–Kutta schemes for hyperbolic systems and kinetic equations in the diffusion limit. SIAM J. Sci. Comput. 35(1), A22–A51 (2013)

    Article  MathSciNet  Google Scholar 

  7. Buet, C., Cordier, S.: An asymptotic preserving scheme for hydrodynamics radiative transfer models. Numer. Math. 108(2), 199–221 (2007)

    Article  MathSciNet  Google Scholar 

  8. Cercignani, C.: The Boltzmann Equation and Its Applications, Applied Mathematical Sciences, vol. 67. Springer, New York (1988)

    Book  Google Scholar 

  9. Cliffe, K.A., Giles, M.B., Scheichl, R., Teckentrup, A.L.: Multilevel Monte Carlo methods and applications to elliptic PDEs with random coefficients. Comput. Vis. Sci. 14(1), 3–15 (2011)

    Article  MathSciNet  Google Scholar 

  10. Crestetto, A., Crouseilles, N., Lemou, M.: A particle micro–macro decomposition based numerical scheme for collisional kinetic equations in the diffusion scaling. Commun. Math. Sci. 16(4), 887–911 (2018)

    Article  MathSciNet  Google Scholar 

  11. Crouseilles, N., Lemou, M.: An asymptotic preserving scheme based on a micro–macro decomposition for collisional Vlasov equations: diffusion and high-field scaling limits. Kinet. Relat. Mod. 4(2), 441–477 (2011)

    Article  MathSciNet  Google Scholar 

  12. Degond, P., Dimarco, G., Pareschi, L.: The moment-guided Monte Carlo method. Int. J. Numer. Methods Fluids 67(2), 189–213 (2011)

    Article  MathSciNet  Google Scholar 

  13. Dimarco, G., Pareschi, L.: Hybrid multiscale methods II. Kinetic equations. Multiscale Model. Simul. 6(4), 1169–1197 (2008)

    Article  MathSciNet  Google Scholar 

  14. Dimarco, G., Pareschi, L.: Fluid solver independent hybrid methods for multiscale kinetic equations. SIAM J. Sci. Comput. 32(2), 603–634 (2010)

    Article  MathSciNet  Google Scholar 

  15. Dimarco, G., Pareschi, L.: High order asymptotic-preserving schemes for the Boltzmann equation. C. R. Math. 350(9–10), 481–486 (2012)

    Article  MathSciNet  Google Scholar 

  16. Dimarco, G., Pareschi, L.: Numerical methods for kinetic equations. Acta Numer. 23, 369–520 (2014)

    Article  MathSciNet  Google Scholar 

  17. Dimarco, G., Pareschi, L., Samaey, G.: Asymptotic-preserving Monte Carlo methods for transport equations in the diffusive limit. SIAM J. Sci. Comput. 40(1), A504–A528 (2018)

    Article  MathSciNet  Google Scholar 

  18. Giles, M.B.: Multilevel Monte Carlo path simulation. Oper. Res. 56(3), 607–617 (2008)

    Article  MathSciNet  Google Scholar 

  19. Giles, M.B.: Multilevel Monte Carlo methods. Acta Numer. 24, 259–328 (2015)

    Article  MathSciNet  Google Scholar 

  20. Goldstein, S.: On diffusion by discontinuous movements, and on the telegraph equation. Q. J. Mech. Appl. Math. 4(2), 129–156 (1951)

    Article  MathSciNet  Google Scholar 

  21. Gosse, L., Toscani, G.: An asymptotic-preserving well-balanced scheme for the hyperbolic heat equations. C. R. Math. 334(4), 337–342 (2002)

    Article  MathSciNet  Google Scholar 

  22. Haji-Ali, A.L., Nobile, F., Tempone, R.: Multi-index Monte Carlo: when sparsity meets sampling. Numer. Math. 132(4), 767–806 (2016)

    Article  MathSciNet  Google Scholar 

  23. Heinrich, S.: Multilevel Monte Carlo methods. In: Lecture Notes in Computer Science (including subseries Lecture Notes in Artificial Intelligence and Lecture Notes in Bioinformatics), vol. 2179, pp. 58–67. Springer (2001)

  24. Hoel, H., Law, K.J.H., Tempone, R.: Multilevel ensemble Kalman filtering. SIAM J. Numer. Anal. 54(3), 1813–1839 (2016)

    Article  MathSciNet  Google Scholar 

  25. Jin, S.: Efficient asymptotic-preserving (AP) schemes for some multiscale kinetic equations. SIAM J. Sci. Comput. 21(2), 441–454 (1999)

    Article  MathSciNet  Google Scholar 

  26. Jin, S., Pareschi, L., Toscani, G.: Diffusive relaxation schemes for multiscale discrete-velocity kinetic equations. SIAM J. Numer. Anal. 35(6), 2405–2439 (1998)

    Article  MathSciNet  Google Scholar 

  27. Jin, S., Pareschi, L., Toscani, G.: Uniformly accurate diffusive relaxation schemes for multiscale transport equations. SIAM J. Numer. Anal. 38(3), 913–936 (2000)

    Article  MathSciNet  Google Scholar 

  28. Klar, A.: An asymptotic-induced scheme for nonstationary transport equations in the diffusive limit. SIAM J. Numer. Anal. 35(3), 1073–1094 (1998)

    Article  MathSciNet  Google Scholar 

  29. Klar, A.: A numerical method for kinetic semiconductor equations in the drift-diffusion limit. SIAM J. Sci. Comput. 20(5), 1696–1712 (1999)

    Article  MathSciNet  Google Scholar 

  30. Lapeyre, B., Pardoux, É., Sentis, R., Craig, A.W., Craig, F.: Introduction to Monte Carlo Methods for Transport and Diffusion Equations, vol. 6. Oxford University Press, Oxford (2003)

    MATH  Google Scholar 

  31. Larsen, E.W., Keller, J.B.: Asymptotic solution of neutron transport problems for small mean free paths. J. Math. Phys. 15(1), 75–81 (1974)

    Article  MathSciNet  Google Scholar 

  32. Lemou, M., Mieussens, L.: A new asymptotic preserving scheme based on micro–macro formulation for linear kinetic equations in the diffusion limit. SIAM J. Sci. Comput. 31(1), 334–368 (2008)

    Article  MathSciNet  Google Scholar 

  33. Løvbak, E., Samaey, G., Vandewalle, S.: A multilevel Monte Carlo asymptotic-preserving particle method for kinetic equations in the diffusion limit. In: L’Ecuyer, P., Tuffin, B. (eds.) Monte Carlo and Quasi-Monte Carlo Methods 2018, pp. 383–402. Springer (2020)

  34. Mortier, B., Baelmans, M., Samaey, G.: Kinetic-diffusion asymptotic-preserving Monte Carlo algorithms for plasma edge neutral simulation. Contrib. Plasma Phys. 60, e201900134 (2019)

    Article  Google Scholar 

  35. Naldi, G., Pareschi, L.: Numerical schemes for hyperbolic systems of conservation laws with stiff diffusive relaxation. SIAM J. Numer. Anal. 37(4), 1246–1270 (2000)

    Article  MathSciNet  Google Scholar 

  36. Pareschi, L.: Hybrid multiscale methods for hyperbolic and kinetic problems. ESAIM Proc. 15, 87–120 (2005)

    Article  MathSciNet  Google Scholar 

  37. Pareschi, L., Caflisch, R.E.: An implicit Monte Carlo method for rarefied gas dynamics. J. Comput. Phys. 154(1), 90–116 (1999)

    Article  MathSciNet  Google Scholar 

  38. Pareschi, L., Russo, G.: An introduction to Monte Carlo method for the Boltzmann equation. ESAIM Proc. 10, 35–75 (2001)

    Article  MathSciNet  Google Scholar 

  39. Pareschi, L., Trazzi, S.: Numerical solution of the Boltzmann equation by time relaxed Monte Carlo (TRMC) methods. Int. J. Numer. Methods Fluids 48(9), 947–983 (2005)

    Article  MathSciNet  Google Scholar 

  40. Pope, S.B.: A Monte Carlo method for the PDF equations of turbulent reactive flow. Combust. Sci. Technol. 25(5–6), 159–174 (1981)

    Article  Google Scholar 

  41. Rousset, M., Samaey, G.: Simulating individual-based models of bacterial chemotaxis with asymptotic variance reduction. Math. Models Methods Appl. Sci. 23(12), 2155–2191 (2011)

    Article  MathSciNet  Google Scholar 

  42. Van Barel, A., Vandewalle, S.: Robust optimization of PDEs with random coefficients using a multilevel Monte Carlo method. SIAM/ASA J. Uncertain. Quantif. 7(1), 174–202 (2019)

    Article  MathSciNet  Google Scholar 

Download references

Acknowledgements

We thank the anonymous reviewer for their substantial suggestions which have significantly improved this article. Emil Løvbak is funded by the Research Foundation - Flanders (FWO) under SB-Fellowship number 1SB1919N. The computational resources and services used in this work were provided by the VSC (Flemish Supercomputer Center), funded by the Research Foundation - Flanders (FWO) and the Flemish Government – department EWI.

Author information

Authors and Affiliations

Authors

Corresponding author

Correspondence to Emil Løvbak.

Ethics declarations

Conflict of interest

The authors declare that they have no conflict of interest.

Additional information

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

A preliminary version of this paper, cited as [33], was published in Monte Carlo and Quasi-Monte Carlo Methods 2018

Supplementary Information

Below is the link to the electronic supplementary material.

Supplementary material 1 (pdf 238 KB)

Appendices

Computing the covariance sums needed to prove Lemma 5

1.1 \(\displaystyle \sum _{m=0}^{M-1} {{\,\mathrm{\mathrm {Cov}}\,}}\left( \Delta T^{n,m}_{p,\ell }, \Delta T^{n}_{p,\ell -1} \right) \)

The covariance between the coarse increment n and a fine sub-increment (nm), given Lemma 4, can be written out as

$$\begin{aligned} {{\,\mathrm{\mathrm {Cov}}\,}}\left( \Delta T^{n,m}_{p,\ell }, \Delta T^{n}_{p,\ell -1} \right) = \frac{\Delta t_{\ell -1}^2 {\tilde{v}}_{\Delta t_\ell } {\tilde{v}}_{\Delta t_{\ell -1}}}{M} {{\,\mathrm{{\mathbb {E}}}\,}}\left[ {\bar{V}}^{n,m}_{p,\ell }{\bar{V}}^{n}_{p,\ell -1}\right] . \end{aligned}$$
(48)

To calculate the expected value in (48), we need to consider the probabilities of coupled collisions taking place in the correlated simulations. If a collision takes place in a given set of M fine increments, the probability that the correlated coarse simulation time step will also simulate a collision is given by

$$\begin{aligned}&P\left( u^n_{p,\ell -1} \ge p_{nc,{\varDelta } t_{\ell -1}}\big | u^{n,\text {max}}_{p,\ell } \ge p_{nc,{\varDelta } t_\ell }\right) \nonumber \\&\quad = P\left( {\left( u^{n,\text {max}}_{p,\ell }\right) }^M \ge p_{nc,{\varDelta } t_{\ell -1}}\big | {\left( u^{n,\text {max}}_{p,\ell }\right) }^M \ge p_{nc,{\varDelta } t_\ell }^M \right) \nonumber \\&\quad = \frac{1-p_{nc,{\varDelta } t_{\ell -1}}}{1-p_{nc,{\varDelta } t_\ell }^M}. \end{aligned}$$
(49)

From Lemma 1, we know that it is not possible for a collision to take place in the coarse simulation, without a collision taking place in the fine simulation. This leaves three possibilities, when considering collision behavior in coarse time step \(n-1\):

  • Both at level \(\ell -1\) and at level \(\ell \), no collision occurred in time step \(n-1\). In this case, time step \(n-1\) will not affect the correlation of the velocities between the simulations. If the velocities were correlated at the beginning of time step \(n-1\), they will still be so at the end of the time step, and vice versa. In this case, we thus need to look at step \(n-2\), and so on, until we reach a past time step that satisfies one of the following two cases.

  • A collision occurred at level \(\ell \) in time step \(n-1\), but not at level \(\ell -1\). In this case, a new \({\bar{V}}^{n-1,m,*}_{p,\ell }\) was drawn from \({{\mathcal {M}}}(v)\), for some m, independently of the value of \({\bar{V}}^{n-1,*}_{p,\ell -1}\). Because all sampled velocities are independent, there is no correlation between \({\bar{V}}^{n,m-1}_{p,\ell }\) and \({\bar{V}}^{n-1}_{p,\ell -1}\), making the expected value of their product zero by (10).

  • A collision occurred both at level \(\ell \) and at level \(\ell -1\) in time step \(n-1\). We know that \({\bar{V}}^{n,0}_{p,\ell } = {\bar{V}}^{n-1,M,*}_{p,\ell } = {\bar{V}}^{n}_{p,\ell -1} = {\bar{V}}^{n-1,*}_{p,\ell -1}\) by (29). For every other fine simulation sub-step \(m>1\) in time step n, there is a probability \(1-p_{nc,{\varDelta } t_\ell }^{m}\) that a collision has taken place in at least one of the steps (n, 0) through (nm). These collisions will not affect the coarse simulation until time step \(n+1\), where \({\bar{V}}^{n,*}_{p,\ell }\) is used to generate a new velocity. So we have that \({\bar{V}}^{n,m}_{p,\ell } = {\bar{V}}^{n}_{p,\ell -1}\) with a probability \(p_{nc,{\varDelta } t_\ell }^{m}\), otherwise the two velocities are uncorrelated.

From the above list of possibilities, we conclude that the random variables \({\bar{V}}^{n,m}_{p,\ell }\) and \({\bar{V}}^{n}_{p,\ell -1}\) are equal if both of the following are true:

  1. 1.

    The last simulation step at level \(\ell \) that underwent a collision is a sub-step of a simulation step \(n^\prime \) at level \(\ell -1\) which also underwent a collision.

  2. 2.

    The coarse simulation step \(n^\prime \) is not the current coarse step n.

Otherwise, \({\bar{V}}^{n,m}_{p,\ell }\) and \({\bar{V}}^{n}_{p,\ell -1}\) are uncorrelated.

Making use of (10) and (49), we have that

$$\begin{aligned} {{\,\mathrm{{\mathbb {E}}}\,}}\left[ {\bar{V}}^{n,m}_{p,\ell }{\bar{V}}^{n}_{p,\ell -1}\right] = \frac{1-p_{nc,{\varDelta } t_{\ell -1}}}{1-p_{nc,{\varDelta } t_\ell }^M} p_{nc,{\varDelta } t_\ell }^{m}. \end{aligned}$$
(50)

Plugging (50) into (48) gives

$$\begin{aligned} \sum _{m=0}^{M-1} {{\,\mathrm{\mathrm {Cov}}\,}}\left( \Delta T^{n,m}_{p,\ell }, \Delta T^{n}_{p,\ell -1} \right)&= \sum _{m=0}^{M-1} \frac{\Delta t_{\ell -1}^2 {\tilde{v}}_{\Delta t_\ell } {\tilde{v}}_{\Delta t_{\ell -1}}}{M} \frac{1-p_{nc,{\varDelta } t_{\ell -1}}}{1-p_{nc,{\varDelta } t_\ell }^M} p_{nc,{\varDelta } t_\ell }^{m} \nonumber \\&= \frac{\Delta t_{\ell -1}^2 {\tilde{v}}_{\Delta t_\ell } {\tilde{v}}_{\Delta t_{\ell -1}}}{M} \frac{p_{c,{\varDelta } t_{\ell -1}}}{p_{c,{\varDelta } t_\ell }} \nonumber \\&=\Delta t_{\ell -1}^2{\tilde{v}}_{\Delta t_{\ell -1}}^2 = {{\,\mathrm{{\mathbb {V}}}\,}}\left[ \Delta T^{n}_{p,\ell -1}\right] , \end{aligned}$$
(51)

where the final equality is a fortuitous coincidence that will shorten some expressions in Sect. 5.

1.2 \(\displaystyle \sum _{m=0}^{M-2} \sum _{m^\prime = m+1}^{M-1}{{\,\mathrm{\mathrm {Cov}}\,}}\left( \Delta T^{n,m}_{p,\ell }, \Delta T^{n,m^\prime }_{p,\ell } \right) \)

As in Sect. A.1, we start by writing out the covariance between subsequent fine sub-increments (nm) and \((n,m^\prime )\) as

$$\begin{aligned} {{\,\mathrm{\mathrm {Cov}}\,}}\left( \Delta T^{n,m}_{p,\ell }, \Delta T^{n,m^\prime }_{p,\ell } \right) = \Delta t_\ell ^2 {\tilde{v}}_{\Delta t_\ell }^2 {{\,\mathrm{{\mathbb {E}}}\,}}\left[ {\bar{V}}^{n,m}_{p,\Delta t_\ell } {\bar{V}}^{n,m^\prime }_{p,\Delta t_\ell } \right] . \end{aligned}$$
(52)

Calculating the sum of the covariance between subsequent fine increments is straightforward. To simplify notation, we introduce \(\Delta m\) as shorthand for \(m^\prime - m\). As the collision probability is constant across time steps, the covariance is given by

$$\begin{aligned} \sum _{m=0}^{M-2} \sum _{m^\prime = m+1}^{M-1} {{\,\mathrm{\mathrm {Cov}}\,}}\left( \Delta T^{n,m}_{p,\ell }, \Delta T^{n,m^\prime }_{p,\ell } \right) = \Delta t_\ell ^2 {\tilde{v}}_{\Delta t_\ell }^2 \sum _{m=0}^{M-2} \sum _{m^\prime = m+1}^{M-1} {{\,\mathrm{{\mathbb {E}}}\,}}\left[ {\bar{V}}^{n,m}_{p,\ell }{\bar{V}}^{n,m^\prime }_{p,\ell }\right] . \end{aligned}$$
(53)

Making use of (10), we get that the r.h.s. of (53) is equal to

$$\begin{aligned} \Delta t_\ell ^2 {\tilde{v}}_{\Delta t_\ell }^2 \sum _{m=0}^{M-2} \sum _{m^\prime = m+1}^{M-1} p_{nc,{\varDelta } t_\ell }^{m^\prime - m} = \Delta t_\ell ^2 {\tilde{v}}_{\Delta t_\ell }^2 \sum _{\mathrel {{\varDelta } m}=1}^{M-1} (M-\mathrel {{\varDelta } m}) p_{nc,{\varDelta } t_\ell }^{\mathrel {{\varDelta } m}} \end{aligned}$$
(54)

By making use of the identity

$$\begin{aligned} \sum _{m=1}^{M-1}(M-m)a^m = \frac{a(a^{M}+M(1-a)-1)}{(1-a)^2} \end{aligned}$$
(55)

and noting that \(\frac{p_{nc,{\varDelta } t_\ell }}{p_{c,{\varDelta } t_\ell }}=\frac{\varepsilon ^2}{\Delta t_\ell }\), we work out (54) as

$$\begin{aligned} \Delta t_\ell ^2 {\tilde{v}}_{\Delta t_\ell }^2 \sum _{\mathrel {{\varDelta } m}=1}^{M-1} (M-\mathrel {{\varDelta } m}) p_{nc,{\varDelta } t_\ell }^{\mathrel {{\varDelta } m}}&= \varepsilon ^2 \Delta t_\ell {\tilde{v}}_{\Delta t_\ell }^2 \frac{p_{nc,{\varDelta } t_\ell }^{M}+Mp_{c,{\varDelta } t_\ell }-1}{p_{c,{\varDelta } t_\ell }} \nonumber \\&= \varepsilon ^2 {\tilde{v}}_{\Delta t_\ell }^2 \left( M\Delta t_\ell - (\varepsilon ^2+\Delta t_\ell )\left( 1 - p_{nc,{\varDelta } t_\ell }^M \right) \right) . \end{aligned}$$
(56)

1.3 \(\displaystyle \sum _{n=0}^{N-2} \sum _{n^\prime =n+1}^{N-1}{{\,\mathrm{\mathrm {Cov}}\,}}\left( \Delta _{T,n} , \Delta _{T,n^\prime } \right) \)

Making use of Lemma 4, we write out the covariance of differences of transport increments for a given n and \(n^\prime \) as

$$\begin{aligned} {{\,\mathrm{\mathrm {Cov}}\,}}\left( \Delta _{T,n},\Delta _{T,n^\prime }\right)&= {{\,\mathrm{{\mathbb {E}}}\,}}\left[ \left( \sum _{m=0}^{M-1} \Delta T^{n,m}_{p,\ell } - \Delta T^n_{p,\ell -1} \right) \left( \sum _{m^\prime =0}^{M-1} \Delta T^{n^\prime ,m^\prime }_{p,\ell } - \Delta T^{n^\prime }_{p,\ell -1} \right) \right] \nonumber \\&= \underbrace{{{\,\mathrm{{\mathbb {E}}}\,}}\left[ \sum _{m=0}^{M-1} \sum _{m^\prime =0}^{M-1} \Delta T^{n,m}_{p,\ell } \Delta T^{n^\prime ,m^\prime }_{p,\ell } \right] }_{(I)} + \underbrace{{{\,\mathrm{{\mathbb {E}}}\,}}\left[ \Delta T^n_{p,\ell -1} \Delta T^{n^\prime }_{p,\ell -1} \right] }_\text {(II)}\nonumber \\&\quad - \underbrace{{{\,\mathrm{{\mathbb {E}}}\,}}\left[ \sum _{m=0}^{M-1} \Delta T^{n,m}_{p,\ell } \Delta T^{n^\prime }_{p,\ell -1} \right] }_\text {(III)} - \underbrace{{{\,\mathrm{{\mathbb {E}}}\,}}\left[ \sum _{m=0}^{M-1} \Delta T^n_{p,\ell -1} \Delta T^{m^\prime ,n^\prime }_{p,\ell } \right] }_\text {(IV)}. \end{aligned}$$
(57)

We now calculate each term on the r.h.s. of (57) separately. We start by writing the expressions in terms of \(\Delta n\), as a shorthand for \(n^\prime - n\). In the remainder of the section we assume \(n^\prime > n\), without loss of generality.

The increment covariance at level \(\ell \) \((\text {I})\) is calculated from a similar starting point to (52):

$$\begin{aligned} {{\,\mathrm{{\mathbb {E}}}\,}}\left[ \sum _{m=0}^{M-1} \sum _{m^\prime =0}^{M-1} \Delta T^{n,m}_{p,\ell } \Delta T^{n^\prime ,m^\prime }_{p,\ell } \right]&= \sum _{m=0}^{M-1} \sum _{m^\prime =0}^{M-1} \Delta t_\ell ^2 {\tilde{v}}_{\Delta t_\ell }^2 {{\,\mathrm{{\mathbb {E}}}\,}}\left[ {\bar{V}}^{n,m}_{p,\ell }{\bar{V}}^{n^\prime ,m^\prime }_{p,\ell }\right] \nonumber \\&= \sum _{m=0}^{M-1} \sum _{m^\prime =0}^{M-1} \Delta t_\ell ^2 {\tilde{v}}_{\Delta t_\ell }^2 p_{nc,{\varDelta } t_\ell }^{M\mathrel {{\varDelta } n}+m^\prime - m}. \end{aligned}$$
(58)

We now substitute the double summation in (58) with a two summations over \(\mathrel {{\varDelta } m}= |m^\prime - m|\):

$$\begin{aligned} \Delta t_\ell ^2 {\tilde{v}}_{\Delta t_\ell }^2&\left( \right. \underbrace{\sum _{\mathrel {{\varDelta } m}=0}^{M-1} (M-\mathrel {{\varDelta } m}) p_{nc,{\varDelta } t_\ell }^{M\mathrel {{\varDelta } n}-\mathrel {{\varDelta } m}}}_{m^\prime \le m} + \underbrace{\sum _{\mathrel {{\varDelta } m}=1}^{M-1} (M-\mathrel {{\varDelta } m}) p_{nc,{\varDelta } t_\ell }^{M\mathrel {{\varDelta } n}+\mathrel {{\varDelta } m}}}_{m^\prime > m} \left. \right) \nonumber \\&= \Delta t_\ell ^2 {\tilde{v}}_{\Delta t_\ell }^2 \left( M + \sum _{\mathrel {{\varDelta } m}=1}^{M-1} (M-\mathrel {{\varDelta } m}) \left( p_{nc,{\varDelta } t_\ell }^{\mathrel {{\varDelta } m}} + p_{nc,{\varDelta } t_\ell }^{-\mathrel {{\varDelta } m}} \right) \right) p_{nc,{\varDelta } t_\ell }^{M\mathrel {{\varDelta } n}}. \end{aligned}$$
(59)

We now again make use of the identity (55), giving that (59) equals

$$\begin{aligned}&\Delta t_\ell ^2 {\tilde{v}}_{\Delta t_\ell }^2 \!\! \left( \!\! M + \frac{p_{nc,{\varDelta } t_\ell }}{(1-p_{nc,{\varDelta } t_\ell })^2}\!\left( \! M\!\left( 2 - p_{nc,{\varDelta } t_\ell }\! -p_{nc,{\varDelta } t_\ell }^{-1} \! \right) \! + p_{nc,{\varDelta } t_\ell }^M \! + p_{nc,{\varDelta } t_\ell }^{-M} \! -2 \right) \!\! \right) p_{nc,{\varDelta } t_\ell }^{M\mathrel {{\varDelta } n}} \nonumber \\&\quad = \varepsilon ^2 {\tilde{v}}^2 \left( p_{nc,{\varDelta } t_\ell }^{1-M} + p_{nc,{\varDelta } t_\ell }^{M+1} - 2 p_{nc,{\varDelta } t_\ell }\right) p_{nc,{\varDelta } t_\ell }^{M\mathrel {{\varDelta } n}}. \end{aligned}$$
(60)

Calculating the covariance of increments at level \(\ell -1\) (II) is also straightforward:

$$\begin{aligned} {{\,\mathrm{{\mathbb {E}}}\,}}\left[ \Delta T^n_{p,\ell -1} \Delta T^{n^\prime }_{p,\ell -1} \right] = \Delta t_{\ell -1}^2{\tilde{v}}_{\Delta t_{\ell -1}}^2 {{\,\mathrm{{\mathbb {E}}}\,}}\left[ {\bar{V}}^n_{p,\ell -1}{\bar{V}}^{n^\prime }_{p,\ell -1}\right] = \Delta t_{\ell -1}^2{\tilde{v}}_{\Delta t_{\ell -1}}^2 p_{nc,{\varDelta } t_{\ell -1}}^{\mathrel {{\varDelta } n}}. \end{aligned}$$
(61)

The expression for the covariance between increments at level \(\ell \) and level \(\ell -1\) at differing time steps n and \(n'\) depends on the relative position of the increments, i.e., whether the increment at level \(\ell \) comes before that at level \(\ell -1\), or not. If the fine increment at level \(\ell \) comes first (III), we need to calculate

$$\begin{aligned} {{\,\mathrm{{\mathbb {E}}}\,}}\left[ \sum _{m=0}^{M-1} \Delta T^{n,m}_{p,\ell } \Delta T^{n^\prime }_{p,\ell -1} \right] = \Delta t_\ell \Delta t_{\ell -1} {\tilde{v}}_{\Delta t_\ell } {\tilde{v}}_{\Delta t_{\ell -1}} \sum _{m=0}^{M-1} {{\,\mathrm{{\mathbb {E}}}\,}}\left[ {\bar{V}}^{n,m}_{p,\ell }{\bar{V}}^{n^\prime }_{p,\ell -1}\right] . \end{aligned}$$
(62)

To calculate the expected value of the r.h.s. of (62), we list the possible simulation behaviors at level \(\ell -1\) in time step n, relative to the fine sub-step m:

  • No collision occurs in the simulation at level \(\ell \). This situation occurs with probability \(p_{nc,{\varDelta } t_\ell }^M.\) In this case, we have already established that no collision occurs at level \(\ell -1\). The probability of the simulations being correlated at the end of time step n is therefore equal to the probability of them being correlated at the start of time step n, given by (49).

  • No collision occurs in sub-steps 0 through \(m-1\), but at least one collision occurs in sub-steps m through \(M-1\). This situation occurs with probability

    $$\begin{aligned} p_{nc,{\varDelta } t_\ell }^{m} \left( 1 - p_{nc,{\varDelta } t_\ell }^{M-m} \right) = p_{nc,{\varDelta } t_\ell }^{m} - p_{nc,{\varDelta } t_\ell }^{M}. \end{aligned}$$

    If the simulation at level \(\ell -1\) also has a collision, then \({\bar{V}}^{n+1}_{p,\ell -1}\) and \({\bar{V}}^{n,m}_{p,\ell }\) are independent by (29). If no collision occurred in the simulation at level \(\ell -1\), which is the case with probability

    $$\begin{aligned} \frac{p_{nc,{\varDelta } t_{\ell -1}}-p_{nc,{\varDelta } t_\ell }^M}{1-p_{nc,{\varDelta } t_\ell }^M}, \end{aligned}$$

    then \({\bar{V}}^n_{p,\ell -1} = {\bar{V}}^{n,m}_{p,\ell }\), if the trajectories were correlated at the beginning of time step n. This is the case with probability (49).

  • At least one collision occurs in sub-steps 0 through \(m-1\), but none occur in sub-steps m through M. This occurs with probability

    $$\begin{aligned} \left( 1 - p_{nc,{\varDelta } t_\ell }^{m} \right) p_{nc,{\varDelta } t_\ell }^{M-m} = p_{nc,{\varDelta } t_\ell }^{M-m} - p_{nc,{\varDelta } t_\ell }^M. \end{aligned}$$

    By (29) we know that \({\bar{V}}^{n,*}_{p,\ell -1} = {\bar{V}}^{n,m}_{p,\ell }\), if a collision also occurs in time step n of the simulation at level \(\ell -1\). This happens with probability (49).

  • Collisions happen both before sub-step m and during or afterwards. In this case, there is no correlation between \({\bar{V}}^{n,*}_{p,\ell -1}\) and \({\bar{V}}^{n,m}_{p,\ell }\).

By adding the non-zero contributions from these four cases, we calculate the probability that \({\bar{V}}^{n,m}_{p,\ell }\) is correlated with \({\bar{V}}^{n,*}_{p,\ell -1}\). To get the probability of the correlation holding until \({\bar{V}}^{n^\prime }_{p,\ell -1}\), we multiply the sum by \(p_{nc,{\varDelta } t_{\ell -1}}^{\mathrel {{\varDelta } n}-1}\). Given the properties (10), we can state that the right hand side sum in (62) is given by

$$\begin{aligned} \sum _{m=0}^{M-1} {{\,\mathrm{{\mathbb {E}}}\,}}&\left[ {\bar{V}}^{n,m}_{p,\ell }{\bar{V}}^{n^\prime }_{p,\ell -1}\right] = \sum _{m=0}^{M-1} \left( p_{nc,{\varDelta } t_\ell }^M + \left( p_{nc,{\varDelta } t_\ell }^{m} \!\! - \! p_{nc,{\varDelta } t_\ell }^{M} \right) \frac{p_{nc,{\varDelta } t_{\ell -1}}-p_{nc,{\varDelta } t_\ell }^M}{1-p_{nc,{\varDelta } t_\ell }^M} \right. \nonumber \\&\quad \left. + p_{nc,{\varDelta } t_\ell }^{M-m} - p_{nc,{\varDelta } t_\ell }^M \right) \frac{1-p_{nc,{\varDelta } t_{\ell -1}}}{1-p_{nc,{\varDelta } t_\ell }^M} p_{nc,{\varDelta } t_{\ell -1}}^{\mathrel {{\varDelta } n}-1}. \end{aligned}$$
(63)

Making us of the fact that \(\frac{p_{c,{\varDelta } t_{\ell -1}}}{p_{nc,{\varDelta } t_{\ell -1}}}=\frac{\Delta t_{\ell -1}}{\varepsilon ^2}\) we can then write the r.h.s. of (63) as

$$\begin{aligned}&\sum _{m=0}^{M-1} \frac{\Delta t_{\ell -1}}{\varepsilon ^2} \! \left( \frac{p_{nc,{\varDelta } t_{\ell -1}}-p_{nc,{\varDelta } t_\ell }^M}{1-p_{nc,{\varDelta } t_\ell }^M} \left( p_{nc,{\varDelta } t_\ell }^{m} \!\! - \! p_{nc,{\varDelta } t_\ell }^{M} \right) + p_{nc,{\varDelta } t_\ell }^{M-m} \right) \!\! \frac{p_{nc,{\varDelta } t_{\ell -1}}^{\mathrel {{\varDelta } n}}}{1-p_{nc,{\varDelta } t_\ell }^M}\nonumber \\&\quad = \frac{\Delta t_{\ell -1}}{\varepsilon ^2} \! \left( \frac{p_{nc,{\varDelta } t_{\ell -1}}-p_{nc,{\varDelta } t_\ell }^M}{1-p_{nc,{\varDelta } t_\ell }^M} \left( \frac{1-p_{nc,{\varDelta } t_\ell }^M}{1-p_{nc,{\varDelta } t_\ell }} - M p_{nc,{\varDelta } t_\ell }^{M} \right) - \frac{p_{nc,{\varDelta } t_\ell }^{M}-1}{p_{nc,{\varDelta } t_\ell }^{-1}-1} \right) \!\! \frac{p_{nc,{\varDelta } t_{\ell -1}}^{\mathrel {{\varDelta } n}}}{1-p_{nc,{\varDelta } t_\ell }^M} \nonumber \\&\quad = \frac{\Delta t_{\ell -1}}{\varepsilon ^2} \left( \frac{ p_{nc,{\varDelta } t_{\ell -1}}-p_{nc,{\varDelta } t_\ell }^M }{ 1-p_{nc,{\varDelta } t_\ell }^M } \left( \frac{1}{p_{c,{\varDelta } t_\ell }} - \frac{Mp_{nc,{\varDelta } t_\ell }^{M}}{ 1-p_{nc,{\varDelta } t_\ell }^M} \right) + \frac{\varepsilon ^2}{\Delta t_\ell } \right) p_{nc,{\varDelta } t_{\ell -1}}^{\mathrel {{\varDelta } n}}. \end{aligned}$$
(64)

To calculate the covariance of coarse time steps preceding fine time steps (IV), fewer calculations are needed. The two time steps are correlated if the trajectories were correlated at the start of coarse simulation time step n and this correlation was not lost due to a collision in the fine simulation between time steps (n, 0) through \((n^\prime ,m^\prime -1)\)

$$\begin{aligned}&{{\,\mathrm{{\mathbb {E}}}\,}}\left[ \sum _{m=0}^{M-1} \Delta T^n_{p,\ell -1} \Delta T^{n^\prime ,m^\prime }_{p,\ell } \right] = \Delta t_\ell \Delta t_{\ell -1} {\tilde{v}}_{\Delta t_\ell } {\tilde{v}}_{\Delta t_{\ell -1}} \sum _{m^\prime =0}^{M-1} {{\,\mathrm{{\mathbb {E}}}\,}}\left[ {\bar{V}}_{p,\ell -1}^n{\bar{V}}_{p,\ell }^{n^\prime ,m^\prime }\right] \nonumber \\&\quad = \Delta t_\ell \Delta t_{\ell -1} {\tilde{v}}_{\Delta t_\ell } {\tilde{v}}_{\Delta t_{\ell -1}} \sum _{m^\prime =0}^{M-1} \frac{1-p_{nc,{\varDelta } t_{\ell -1}}}{1-p_{nc,{\varDelta } t_\ell }^M} p_{nc,{\varDelta } t_\ell }^{M\mathrel {{\varDelta } n}+ m^\prime } \nonumber \\&\quad = \Delta t_\ell \Delta t_{\ell -1} {\tilde{v}}_{\Delta t_\ell } {\tilde{v}}_{\Delta t_{\ell -1}} \frac{1-p_{nc,{\varDelta } t_{\ell -1}}}{1-p_{nc,{\varDelta } t_\ell }^M} \frac{1-p_{nc,{\varDelta } t_\ell }^M}{1-p_{nc,{\varDelta } t_\ell }} p_{nc,{\varDelta } t_\ell }^{M\mathrel {{\varDelta } n}} \nonumber \\&\quad = \Delta t_{\ell -1}^2 {\tilde{v}}_{\Delta t_{\ell -1}}^2 p_{nc,{\varDelta } t_\ell }^{M\mathrel {{\varDelta } n}}. \end{aligned}$$
(65)

Making use of expressions (60) through (65), we get

$$\begin{aligned}&\sum _{n=0}^{N-1} \sum _{n^\prime =n}^N {{\,\mathrm{\mathrm {Cov}}\,}}(\Delta _{T,n},\Delta _{T,n^\prime }) \nonumber \\&\quad = \sum _{n=0}^{N-1} \sum _{n^\prime =n}^N \sigma _1p_{nc,{\varDelta } t_\ell }^{M({n - n^\prime })} + \sigma _2p_{nc,{\varDelta } t_{\ell -1}}^{{n - n^\prime }}\nonumber \\&\quad = \sum _{\mathrel {{\varDelta } n}=1}^{N-1} (N-\mathrel {{\varDelta } n}) \left( \sigma _1p_{nc,{\varDelta } t_\ell }^{M\mathrel {{\varDelta } n}} + \sigma _2p_{nc,{\varDelta } t_{\ell -1}}^{\mathrel {{\varDelta } n}} \right) \nonumber \\&\quad = \sigma _1p_{nc,{\varDelta } t_\ell }^M \!\! \frac{p_{nc,{\varDelta } t_\ell }^{MN} \!\! - \! Np_{nc,{\varDelta } t_\ell }^M \!\! + \! N \! - \! 1}{\left( 1 - p_{nc,{\varDelta } t_\ell }^M \right) ^2} + \sigma _2p_{nc,{\varDelta } t_{\ell -1}}\!\! \frac{p_{nc,{\varDelta } t_{\ell -1}}^{N} \!\!\! - \! Np_{nc,{\varDelta } t_{\ell -1}}\!\! + \! N \! - \! 1}{\left( 1-p_{nc,{\varDelta } t_{\ell -1}}\right) ^2}, \end{aligned}$$
(66)

\( \text {with} \qquad \sigma _1= \varepsilon ^2 {\tilde{v}}^2 \left( p_{nc,{\varDelta } t_\ell }^{1-M} + p_{nc,{\varDelta } t_\ell }^{M+1} - 2 p_{nc,{\varDelta } t_\ell }\right) - \Delta t_{\ell -1}^2 {\tilde{v}}_{\Delta t_{\ell -1}}^2 \qquad \text {and} \)

$$\begin{aligned} \sigma _2= \Delta t_{\ell -1}^2{\tilde{v}}_{\Delta t_{\ell -1}}^2 \!\! \left( \! 1 \! - \! \frac{\Delta t_\ell {\tilde{v}}_{\Delta t_\ell }}{\varepsilon ^2 {\tilde{v}}_{\Delta t_{\ell -1}}} \!\! \left( \frac{ p_{nc,{\varDelta } t_{\ell -1}}-p_{nc,{\varDelta } t_\ell }^M }{ 1-p_{nc,{\varDelta } t_\ell }^M } \! \left( \! \frac{1}{p_{c,{\varDelta } t_\ell }} - \frac{Mp_{nc,{\varDelta } t_\ell }^{M}}{ 1-p_{nc,{\varDelta } t_\ell }^M} \right) \! + \frac{\varepsilon ^2}{\Delta t_\ell } \right) \!\! \right) \!\! . \end{aligned}$$

Theoretical background to Section 6.2

Based on Section 2.6 in [19] it can be shown that it is beneficial to leave out level 0 if

$$\begin{aligned} \sqrt{C_0 {{\,\mathrm{{\mathbb {V}}}\,}}\left[ F_0 \right] } + \sqrt{\left( C_0 + C_1 \right] ) {{\,\mathrm{{\mathbb {V}}}\,}}\left[ F_1-F_0 \right] } - \sqrt{C_1 {{\,\mathrm{{\mathbb {V}}}\,}}\left[ F_1 \right] } > 0, \end{aligned}$$
(67)

with \(C_\ell \) the cost of a sample at level \(\ell \). Our analytical results from Sect. 6.1 can only strictly be applied in (67) if the quantity of interest is the expected particle position or velocity, as applying Lipschitz constants gives independent upper bounds for the positive and negative terms, hence we will also back up our claims with numerical experiments. Here, we only consider the particle position variance, as \({{\,\mathrm{{\mathbb {V}}}\,}}\left[ X_\ell -X_{\ell -1}\right] \gg {{\,\mathrm{{\mathbb {V}}}\,}}\left[ V_\ell -V_{\ell -1}\right] \) and \({{\,\mathrm{{\mathbb {V}}}\,}}\left[ X_\ell \right] \gg {{\,\mathrm{{\mathbb {V}}}\,}}\left[ V_\ell \right] \), for large \(\Delta t_\ell \).

We first fix \(\Delta t_1 = \varepsilon ^2\) and solve (67) for positive real values of \(M>1\) for which the inequality becomes an equality, taking different values of \(\varepsilon \) and \(t^*\). In Fig. 10 we show these unique, numerically computed M-values. It is clear from the range of the color bar that this M-value consistently exists and lies within the interval [6, 13], meaning that any sign change in the r.h.s. of (67) lies in this range. Next, we evaluate the left hand side of (67) for \(M=6\) and \(M=13\) and plot the results in Fig. 11. From this figure, it is clear that (67) holds for small values of M, but no longer holds once the threshold value from Fig. 10 is passed, meaning that (67) consistently holds once M is sufficiently large. We can thus conclude that it makes sense, from a theoretical point of view, to add an extra coarse level, if it is sufficiently coarse in comparison with the level with \(\Delta t = \varepsilon ^2\).

Fig. 10
figure 10

M-value causing a sign change in the left hand side of (67)

Fig. 11
figure 11

Left hand side of (67)

Rights and permissions

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Løvbak, E., Samaey, G. & Vandewalle, S. A multilevel Monte Carlo method for asymptotic-preserving particle schemes in the diffusive limit. Numer. Math. 148, 141–186 (2021). https://doi.org/10.1007/s00211-021-01201-y

Download citation

  • Received:

  • Revised:

  • Accepted:

  • Published:

  • Issue Date:

  • DOI: https://doi.org/10.1007/s00211-021-01201-y

Mathematics Subject Classification

Navigation